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ABSTRACT 


This  thesis  investigates  the  behavior  of  Graphite -Epoxy 
composites  subjected  to  fires  as  may  occur  on  the  decks  of 
naval  aircraft  carriers.  The  analytical  model  consists  of  two 
parts;  one  for  the  determination  of  the  temperature  field 
within  the  composite  due  to  a  fire  and  the  other  for 
determining  the  stresses  within  the  composite  due  to  the 
temperature  field.  Both  problems  are  provided  one -dimensional 
finite  element  models.  Appropriate  failure  criteria  are 
incorporated  to  predict  the  survivability  of  composites  in 
various  fire  environments.  Parametric  studies  were  performed 
with  the  results  presented  in  both  graphical  and  tabular  form. 
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I .  INTRODUCTION 


A.  BACKGROUND 

Composite  materials,  due  to  their  high  strength  to  weight 
ratio,  are  becoming  increasingly  common  as  construction 
materials  in  today's  modern  combat  aircraft.  Through  use  of 
these  materials  the  weight  to  strength  ratio  of  the  airframe 
has  been  greatly  reduced  allowing  for  an  increase  in  the 
respective  payload  per  unit  weight  of  the  structure. 

Several  different  Navy  fixed  wing  aircraft  such  as  the  F- 
18,  F-16,  F-15  and  F-14  all  use  composite  materials  in  their 
construction.  Since  these  aircraft  are  used  on  board  naval 
aircraft  carriers,  it  is  quite  conceivable  that  these 
structures  could  be  exposed  to  flight  deck  fires  of  varying 
intensities  during  their  service  life.  It  is  important  for 
the  Navy  to  assess  the  survivability  of  composite  aircraft 
exposed  to  fires.  To  this  end,  experimental  and  analytical 
investigations  have  been  undertaken. 

In  1979  testing  was  performed  at  the  Naval  Weapons  Center 
(China  Lake)  to  determine  the  relative  survivability  of 
Graphite -Epoxy  Composite  materials  and  7075-T6  Aliiminum  Alloy 
when  exposed  to  flight  deck  pool  fires  [Ref.  l] .  These  tests 
indicated  the  susceptibility  to  burn  through  for  Graphite - 
Epoxy  samples  was  approximately  the  same  as  for  the  comparable 
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7075 -T6  Aluminum  Alloy  scimples  tested.  However  unlike  the 
aluminum,  in  some  cases  the  composite  material  continued  to 
smolder  and/or  burn  after  application  of  extinguishing  agents. 
In  an  attempt  to  model  this  smoldering/burning  of  Graphite- 
Epoxy  Composites,  Vatikiotis  in  1980  [Ref.  2]  developed  a 
computer  model  which  he  later  refined  in  1982  [Ref.  3] .  This 
finite  element  program  evaluates  fiber  diameters,  fiber 
temperatures,  air  temperatures  and  oxygen  levels  throughout 
Graphite -Epoxy  materials  exposed  to  fire.  The  evaluations  are 
with  respect  to  time  and  location. 

Graphite -Epoxy  material  consists  of  high  strength  graphite 
fibers  surrounded  by  an  epoxy  matrix.  When  exposed  to  a  fire 
the  epoxy  matrix  material  burns  off  rather  quickly  at  400- 
500°F  leaving  only  the  graphite  fibers  supported  by  residual 
combustion  products.  These  remaining  high  strength  fibers 
develop  thermally  induced  stresses  due  to  the  varying 
temperature  profile  across  the  material.  If  high  enough 
stresses  are  developed,  fibers  may  fail  rendering  the 
remaining  material  of  no  structural  use.  However,  if  a 
considerable  amount  of  fibers  do  not  fail  the  surviving 
material,  although  considerably  weakened  in  the  compression 
mode,  could  continue  to  be  used  in  the  tensile  mode. 

B.  OBJECTIVE  OF  THIS  THESIS 

This  thesis  investigated  thermally  induced  stress  failure 
of  Graphite -Epoxy  material  exposed  to  varying  intensity  fires. 
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The  analysis  was  carried  out  in  two  parts.  In  part  one,  the 
computer  model  developed  by  Vatikiotis  in  1982  was  used  to 
model  fires  of  different  intensities.  This  program  evaluated 
the  temperature  profiles  across  a  porous  medium  consisting  of 
graphite  fibers  as  time  progressed.  These  results  were 
further  utilized  in  the  second  part  of  the  analysis.  Part  two 
of  the  analysis  consisted  of  developing  a  stress  model  and 
finite  element  computer  program.  This  stress  code  evaluated 
the  thermally  induced  stresses  of  fibers  based  on  the 
temperature  profiles  determined  from  the  combustion  code.  The 
fiber  stresses  were  calculated  at  each  time  step;  failure 
criteria  applied  and  failed  fibers  removed  from  the  stress 
analysis.  Through  this  procedure,  a  quantitative  evaluation 
of  Graphite -Epoxy  material  failure  with  respect  to  time  was 
obtained. 

C.  THESIS  OUTLINE 

Chapter  II  presents  the  physical  and  mathematical  models. 

Chapter  III  develops  the  failure  criteria  used  in 
determination  of  individual  fiber  failure. 

Chapter  IV  contains  a  description  of  the  stress  analysis 
program  and  it's  interaction  with  the  combustion  program 
developed  by  Vatikiotis. 

Chapter  V  presents  case  studies  with  results  presented  in 
graphical  form. 

Chapter  VI  contains  the  conclusions. 
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II.  PROBLEM  FORMDLRTION 


A.  PHYSICAL  MODEL 

The  physical  model  assiimed  for  this  study  is  that  of  a 
Graphite -Epoxy  composite  plate  exposed  to  fires  of  varying 
intensity  as  shown  in  Figure  1.  This  plate  of  composite 
material  could  have  several  different  nominal  thicknesses 
depending  on  the  nximber  of  laminate  layers  or  ply  assumed. 
For  this  study,  nominal  thicknesses  of  0.25,  0.5,  0.75  and  1.0 
inch  are  assumed,  corresponding  to  approximately  35,  70,  105 
and  140  ply  composites  respectively. 
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As  the  composite  plate  heats  up  it  will  lose  its  epoxy 
matrix  material  between  400  and  500°F.  The  matrix  loss 
process  is  too  complex  to  model  exactly.  However,  for  epoxy 
loss  to  occur,  the  epoxy  matrix  must  undergo  phase  changes 
from  solid  to  combustion  products.  These  phase  changes  occur 
at  a  constant  temperature  throughout  the  material. 

This  study  assumes  that  no  failure  of  graphite  fibers 
occurs  prior  to  the  burn-off  of  the  epoxy  matrix  material.  The 
point  at  which  the  epoxy  matrix  has  been  reduced  to  residue 
will  be  chosen  as  the  starting  time  for  all  case  studies. 
Because  of  the  phase  changes  which  occur  during  the  matrix 
burn  off  process  described  above,  the  initial  temperature 
profile  across  the  remaining  Graphite  fiber  medium  is  assumed 
to  be  a  constant  400°F  at  the  starting  time. 

Air  flows  across  the  upper  surface  of  the  composite  plate. 
This  surface  air  flow  produces  a  differential  pressure  across 
the  thickness  of  the  plate.  At  the  starting  time  for  these 
studies,  since  the  matrix  material  has  been  removed,  this 
differential  pressure  will  induce  air  flow  perpendicular 
through  the  remaining  composite  plate. 

B.  VATIKIOTIS  COMBUSTION  MODEL 

In  order  to  obtain  the  temperature  field  required  for 
stress  analysis  the  one -dimensional  finite  element  program 
developed  by  Vatikiotis  was  used.  As  seen  in  Figure  2  his 
combustion  model  assumes  the  composite  material  has  been 
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exposed  to  heat  sufficient  to  remove  the  epoxy  matrix  leaving 
only  the  graphite  fibers  surrounded  by  air  filled  voids.  The 
heat  supplied  by  the  fire  is  implemented  through  application 
of  a  constant  heat  flux  boundary  condition  at  the  lower,  X=0, 
surface.  Differential  pressure  across  the  plate  thickness 
forces  air  flow  through  the  material.  Prior  to  fiber 
combustion  the  incoming  air  is  considered  to  be  of  standard 
composition  containing  no  fire  combustion  products.  After 
combustion  begins  carbon  dioxide  and  carbon  monoxide  also  flow 
through  the  material . 


Figure  2:  The  Vatikiotis  combustion  model. 

In  his  presentation  Vatikiotis  divides  the  material  across 
its  thickness  into  30  equal  thickness  elements  having  a  total 
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of  31  nodal  points.  He  further  develops  a  system  of  four 
transient  one -dimensional  equations  consisting  of  an  energy 
balance  on  the  graphite  material,  an  energy  balance  on  the 
internal  air,  conservation  of  species  on  oxygen  and  a 
combination  of  Darcy's  law  and  continuity  equation  for 
internal  air  flow  velocities. 

The  two  energy  balances  considered  a  differential  volume 
of  material  with  the  following  convention: 

Heat  into  ^  Heat  ^  Heat  out  of  ^  Increased 
dV  Generation  dV  Internal  Energy 

The  energy  balance  for  the  graphite  fibers  accounted  for 
heat  transfer  into  and  out  of  the  differential  volume  by  three 
modes;  conduction,  convection  and  radiation,  using  first  term 
Taylor  series  approximations.  The  heat  generation  term, 
resulting  from  burning  of  the  graphite  fibers,  was  modeled 
using  an  expression  of  Arrhenius  type. 

The  energy  balance  for  the  internal  air  considers  only 
conduction  and  convection  at  the  differential  voliome  surfaces 
along  with  energy  transport  due  to  the  flow  of  air  through  the 
volume.  Also  this  equation  does  not  contain  a  heat  generation 
term. 

The  oxygen  mass  balance  started  from  the  basic  form  of: 

Oxygen  into  _  Oxygen  out  of  ^  Oxygen  ^  Oxygen 

dV  dV  Consumption  Accumulation 

The  oxygen  into  and  out  of  the  differential  volume 
considers  both  molecular  diffusion  and  convective  transport 
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across  the  differential  volume  surfaces,  again  using  a  first 
term  Taylor  series  approximations.  The  oxygen  consumption  was 
based  on  the  graphite  fibers  (carbon)  burning  with  the  oxygen 
present,  in  a  appropriately  specified  stoichiometric  ratio. 
This  reaction  is  assumed  to  produce  carbon  dioxide  and  carbon 
monoxide  only.  A  relationship  for  the  consumption  of  oxygen 
was  developed  through  combination  of  the  stoichiometric  ratio 
with  the  Arrhenius  reaction  rate  equation. 

The  final  equation  of  this  system  is  the  equation  for 
internal  air  flow  velocities.  This  equation  was  developed 
using  a  combination  of  Darcy's  law  for  flow  in  a  porus  medium 
and  the  continuity  equation. 

The  above  four  equations  form  a  set  of  coupled,  nonlinear, 
transient,  partial  differential  equations.  These  equations 
are  solved  for  graphite  temperature,  air  temperature,  oxygen 
concentration,  and  internal  air  flow  velocities  with  respect 
to  time  and  position.  Incorporated  in  this  mathematical  model 
are  the  variations  of  temperature  dependent  properties  namely, 
graphite  thermal  conductivity,  fluid  (air)  viscosity,  and  air 
density,  along  with  parcuneters  derived  from  them.  These 
properties  and  parameters  are  recalculated  at  each  time  step 
of  the  transient  solution.  The  Darcy's  law- continuity 
equation  was  solved  using  a  shooting  method  and  the  remaining 
three  equations  were  solved  using  a  Galerkin  formulation  of 
the  finite  element  method.  For  completeness  the  above 
equations  in  their  final  form  along  with  their  associated 
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boundary  conditions  are  contained  in  Appendix  A.  For  more  in 
depth  review  the  reader  is  referred  to  Vatikiotis'  original 
works . 

The  main  observations  of  Vatikiotis'  earlier  work  are  very 
obvious  when  his  computer  code  is  used.  The  differential 
pressure  across  the  porous  medium  causes  a  flow  of  air  to  flow 
through  medium.  This  induced  air  flow  provides  for  both 
convective  heat  transfer  of  the  fibers  and  a  ready  supply  of 
oxygen  for  the  combustion  process.  These  two  effects  of  the 
air  flow  are  conflicting  effects.  First,  the  supply  of  oxygen 
increases  the  rate  of  combustion  which  increases  the  fiber 
temperature.  Second,  the  convective  heat  transfer  away  from 
the  hot  fibers  reduces  the  fiber  temperature.  A  system 
comprised  of  a  plate  and  a  fire  is  defined  by  (a)  the  fire 
intensity  and  its  time  duration  and  {b)  the  plate  subjected  to 
the  fire.  A  plate  with  known  initial  temperatures  and  oxygen 
concentrations,  thickness,  porosity,  and  pressure  differential 
is  subjected  to  a  fire  of  specified  intensity  for  a  given  time 
period.  When  the  fire  is  removed  the  plate  may  either 
continue  to  burn  (plate  combustion)  or  the  plate  may  return  to 
ambient  conditions  (plate  extinction) .  Which  process  occurs 
depends  on  the  state  of  the  plate  when  the  fire  has  been 
removed.  If  the  temperature  and  air  velocity  fields  exceed  a 
critical  state  then  the  convective  heat  transfer  from  the 
fibers  to  the  air  is  overtaken  by  the  heat  generated  by  the 
fiber  combustion  and  the  material  will  continue  to  burn.  If 
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the  critical  state  was  not  achieved,  the  convective  heat 
transfer  will  predominate  resulting  in  cooling  and  extinction. 

C.  STRESS  ANALYSIS  MODEL 

The  combustion  model  described  previously  determines  the 
temperature  distribution  across  the  thickness  of  the  medium. 
The  stress  analysis  model  uses  the  temperature  profiles  from 
the  combustion  code  to  determine  the  thermally  induced 
stresses  in  the  fiber  composite  plate. 

As  seen  in  Figure  3  the  stress  model  was  developed 
assuming  the  graphite  fibers  were  cylindrical  rods  running 
parallel  to  the  plate  surfaces.  These  fiber  rods  are  modeled 
as  rigidly  fixed  at  one  end  and  on  the  opposite  end  are 
modeled  to  move  in  unison  with  the  other  fibers.  The 
expansion  or  contraction  of  fibers  due  to  varying  temperatures 
will  cause  forces  to  develop  in  the  other  fibers  (See  Figure 
4) .  These  resulting  forces  could  cause  excessive  stresses  to 
develop  possibly  causing  failure.  The  final  solution  to  this 
statically  indeterminate  problem  required  use  of  finite 
element  techniques. 
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THE  EFFECT  OF  HEATING 


PRIOR  TO  HEATING 


Plgur«  4:  The  effect  of  heating  graphite  rods. 
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INAL 


As  seen  in  Figure  5  each  element  (fibers  per  ply)  was 
modeled  as  a  linear  elastic  spring  with  forces  applied  at  each 
end  (nodal  point).  The  finite  element  "direct  approach"  [Ref. 
7]  was  applied  to  each  individual  element  to  obtain  the 
element  stiffness  matrix  and  force  vector  as  follows.  All 
forces,  displacements,  stresses,  strains  and  element 
characteristics  are  in  the  axial  direction  only. 

[— WWf-r 

I — ►Ui  ’ — ►Uz 

Figure  5:  The  linear  spring  element  of  the  finite  element 
direct  approach. 

First  referring  to  Figure  5,  using  summation  of  forces: 

5]f=Fi  +  F2  =  0  (2.1) 

Fi  =  -F2  (2.1a) 

where : 

Fj  =  force  applied  at  point  1 
Fj  =  force  applied  at  point  2 
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Assuming  a  linear  spring  model  yields: 


F2=k(U2-u^)  (2.2) 

where : 

=  displacement  of  point  1 
U2  =  displacement  of  point  2 
k  =  stiffness  coefficient. 

From  [Ref.  4]  the  stiffness  coefficient  for  a  cylindrical 
rod  with  axial  displacement  is: 


k  = 


AE 

L 


(2.3) 


where : 

A  *  cross  sectional  area  of  the  bar  (fibers  per  ply) 
E  =  elastic  modulus  of  fiber  (axial) 

L  »  length  of  fiber. 

Combining  equations  (2.1a)  and  (2.2)  yields  the  element 
matrix  equations  of: 


■  k 

-k 

^  ' 

K) 

-k 

k 

n 

N 

n 

for  fibers,  where  n  is  the  ply  number.  The  forces  and  F2 
are  equal  in  magnitude  and  come  from  the  thermal  expansion  of 
the  material  as  follows. 
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The  thermal  strain  of  the  material  is  defined  as: 

€,^  =  a(r-rj  (2.5) 

where : 

=  thermal  strain 

a  =  thermal  coefficient  of  expansion 
T  =  the  temperature  of  the  material 
T„  =  the  temperature  at  which  the  strain  is  zero. 
Multiplication  of  equation  (2.5)  by  Youngs  modulus  E  yields 
the  stress  c, 

a=E€=Ea(T-Tj  (2.6) 

which  results  if  free  expansion  or  contraction  is  prevented. 
Multiplying  equation  (2.6)  by  the  cross  sectional  area  of  the 
element  yields  the  element  forces; 

-  (2.7) 

at  the  ends  of  the  elements  where  n  is  the  ply  number. 
Dropping  the  subscripts  1  and  2  on  the  forces  and  combining 
equation  (2,7)  with  (2.4)  and  (2.3)  gives  the  following  matrix 
equations  for  each  element: 

k  -k 
-k  k 

where : 

k  =>  AE/L  based  on  fiber  properties  and  dimensions 
Uj^  &  U2  =  displacements  of  points  1  &  2  for  fiber  n 
F  =  AEa(T-T„)  . 


. ■  0. 


(2.8) 
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Using  the  numbering  system  of  Figure  6,  the  global 
stiffness  matrix  and  global  force  vector  are  assembled  as: 

(2.9) 


'^1 

-■^1 

0 

0 

0 

0 

0 

0 

0 

0 
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• 

>1 
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0 

0 
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• 
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where : 

kjj  =  a  component  in  the  global  stiffness  matrix 

U„  =  a  component  in  the  global  displacement  vector 

Fjj  =  a  component  in  the  global  force  vector. 

The  boundary  conditions  imposed  on  the  above  system  of 
equations  are  as  follows.  For  the  rigidly  fixed  nodal  points, 

U^  =  U,=U^  =  U^= _ =  (^2;,-l  =  0  (2.10) 

and  for  the  nodal  points  that  move  in  unison, 

C7^  =  y,  =  a,  =  U3= - =  C72„  .  (2.11) 

The  system  of  equations  of  equation  (2.9)  with  the 
boundary  conditions  of  (2.10)  and  (2.11)  applied  were  solved 
for  nodal  point  displacements  using  IMSL  math  library 
subroutine  DLSARG. 
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solution. 

The  stresses  of  the  fibers  per  ply  were  calculated  using 
the  displacements  obtained  from  the  above  system  of  equations 
as  follows  (see  figure  4) : 

(2.12) 

where 

«  strain  from  actual  motion  of  nodal  points- (Uj-Uj^) /L 
-  strain  causing  fiber  stress 
“  fiber  strain  for  free  expans ion-a (T-T.)  . 
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Substitution  of  appropriate  expressions  for  and  gives 
an  expression  for  mechanical  strain  of: 


.  (2-13) 

Finally  multiplying  equation  (2.13)  by  Young's  modulus  yields: 

o  =  =  -a  (r-r.)|  (2.14) 

where  a  is  the  stress  developed  in  the  fiber. 

Using  the  preceding  equations  a  computer  stress  program 
was  developed  utilizing  the  Vatikiotis  combustion  code  for 
fiber  temperature  profiles.  The  stress  program  determined 
fiber  stresses  throughout  the  composite  material  as  a  function 
of  time  and  location.  The  stress  program  also  determined 
failed  fibers  basod  on  failure  criteria  of  Chapter  III. 
Specifics  of  the  stress  program  are  contained  in  Chapter  IV. 
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III.  FAILURE  CRITERIA 


A.  INTRODUCTION 

The  stress  analysis  code  determines  fiber  stresses 
throughout  the  composite  plate  and  then  determines  if  failure 
of  any  of  the  fibers  has  occurred.  In  this  chapter  possible 
fiber  failure  criteria  are  considered.  The  four  failure  modes 
considered  were  tensile  strength  failure,  compressive  strength 
failure,  fiber  buckling,  and  fiber  consumption.  Of  the  four 
modes  considered  only  three  were  incorporated  in  the  stress 
analysis  code. 

B.  TENSILE  FAILURE 

The  first  criteria  considered  was  failure  due  to  exceeding 
the  fibers  tensile  strength  of  300  ksi.  When  subjected  to 
heat  graphite  fibers  tend  to  contract,  which  is  reflected  in 
their  thermal  expansion  coefficient  of  -0 . 55X10 
Initially,  the  fiber  closest  to  the  heat  source  will  increase 
in  temperature  then  over  time  the  remaining  fibers  will  heat. 
The  plate  fibers  will  continue  to  rise  in  temperature  as  long 
as  the  heat  source  is  applied.  Since  thermal  strain  is 
proportional  to  temperature,  the  magnitude  of  fiber 
contraction  increases  with  fiber  temperature.  Due  to  the 
fibers  closest  to  the  fire  being  hotter  than  those  further 
from  the  fire,  fibers  closest  to  the  fire  would,  if  free  to  do 
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so,  contract  more  than  those  farther  away.  This  variation  in 
tendency  to  contract  causes  induced  stresses  since  all  fibers 
are  restrained  to  contract  uniformly  as  seen  in  Figure  4 . 


In  order  to  estimate  the  temperature  required  to  reach 
tensile  strength,  a  limiting  case  of  Equation  (2.14)  was  used. 
Since  the  fiber  closest  to  the  heat  source  is  the  first  to 
heat,  its  temperature  begins  to  rise  before  any  significant 
change  in  the  other  fiber  temperatures.  As  a  worst  case,  this 
fiber  could  be  heated  to  a  very  high  temperature  while  all  the 
remaining  fibers  remained  very  close  to  their  initial 
temperature.  If  the  number  of  fibers  were  relatively  high, 
this  situation  could  be  modeled  as  one  heated  fiber  with  fixed 
ends  (i.e.  not  allowed  to  move  axially)  .  Using  this  model 
equation  (2.14)  becomes: 

Sc=-Ea(T-Tj  (3.1) 

or,  by  rearranging,  the  minimizm  temperature  required  for  this 
fiber  to  reach  the  tensile  strength  is; 


Sr 

T  ■  = _ —  +  T 

-‘min 


(3.2) 


where : 

S^.  =  tensile  strength  of  the  fiber  (300  ksi) 

'^min  “  minimum  temperature  for  tensile  failure 
E  =  elastic  modulus  (31.0X10®  psi) 

Tjo  =  reference  temperature  (assumed  to  be  60°F) 
a  =  thermal  expansion  coefficient  (-0.55X10'®/°F)  . 
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Solving  the  above  equation  for  the  minimum  fiber 
temperature  gives  a  value  of  approximately  17,600°F.  This 
temperature  is  unachievable  since  graphites  melting  point  is 
approximately  6000°F  and,  moreover,  the  combustion  code  is 
unable  track  temperatures  over  approximately  2600°F. 

Due  to  the  above  calculation,  the  probability  of  any 
fibers  failing  in  the  tensile  mode  is  highly  unlikely. 
However,  a  quick  check  for  tensile  failure  was  incorporated 
into  the  computer  program  as  verification  of  the  above 
calculations.  Tensile  failure  was  not  observed  during  any  of 
the  case  studies. 

C.  COMPRESSIVE  STRENGTH  FAILURE 

The  next  criteria  considered  was  purely  axial  compression 
of  the  fibers.  For  the  graphite  fibers  under  consideration 
their  compressive  strength  is  260  ksi.  However,  this  high 
value  is  quite  deceiving.  These  fibers  have  a  very  small 
cross  section  and  relatively  long  length,  causing  them  to  be 
susceptible  to  buckling  failure  when  in  compression.  Stresses 
at  buckling  are  much  lower  than  the  compressive  fracture 
strength.  For  this  reason  compressive  strength  of  the  fibers 
was  not  used  as  a  failure  criteria. 

D.  FIBER  BUCKLING 

The  third  failure  mode  considered  was  fiber  buckling.  The 
graphite  fibers  have  a  very  high  slenderness  ratio  due  to 
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their  small  cross  sectional  diameter  and  relatively  long 
length.  These  fibers  could  be  susceptible  to  Euler  column 
buckling.  Evaluation  of  this  mode  of  failure  requires  a 
"typical"  length  for  the  unsupported  fibers.  This  fiber 
length  was  assumed  to  be  10  inches,  corresponding  to  a 
reasonable  sized  hole  burned  through  the  material.  These 
fibers  were  further  assximed  to  have  fixed  end  conditions. 

The  development  of  Euler  Buckling  formula  is  available  in 
several  texts  and  will  not  be  redeveloped  here.  From  [Ref.  5] 
the  formula  for  Euler  column  buckling  assuming  fixed  end 
conditions  with  constant  circular  cross  sectional  area  is: 


a 


CZ 


4L2 


(3.3) 


where : 

=  critical  stress  for  column  buckling 
E  =  elastic  modulus  of  the  fibers 
D  =  the  Diameter  of  the  fibers 
L  =  the  length  of  the  unsupported  fibers. 

Using  the  above  formula,  the  calculated  value  of  cr^j.  is 
approximately  19  psi.  This  value  is  very  small  and  quite 
conservative  since  the  model  assumes  no  lateral  support. 

As  the  epoxy  matrix  burns,  some  of  the  combustion  products 
will  stay  behind  forming  a  new  matrix  material.  This  new 
matrix  structure  consisting  of  residual  combustion  products, 
although  weaker  than  the  original  epoxy  matrix,  renders 
lateral  support  to  the  fibers.  This  additional  support 
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significantly  increases  the  critical  buckling  strength  of  the 
fibers . 

Rosen  [Ref.  6]  ,  in  his  work  on  critical  buckling  stress  of 
composite  materials,  develops  two  analytical  models  for 
buckling  of  composites,  namely;  the  extensional  mode  and  the 
shear  mode.  Of  these  two  modes,  the  mode  most  likely  for 
fiber  volume  fractions  of  approximately  20%  is  the  shear  mode. 
The  shear  mode  formula  is : 


'fez' 


V,{l-V,) 


(3.4) 


where : 

=  fiber  critical  buckling  stress 
Vf  =  fiber  volume  fraction 
G„  *  shear  modulus  of  the  matrix  material . 

A  search  of  literature  provided  little  information  on  the 
properties  of  epoxy  combustion  products  for  use  as  a  matrix 
material.  However,  after  some  consideration,  this  problem  was 
overcome  by  assuming  that  the  effect  of  the  combustion 
products,  as  a  matrix  material,  could  be  modeled  as  a  small 
fraction  (F)  of  the  original  graphite- epoxy  composite  buckling 
critical  stress  given  by  equation  (3.4)  .  This  fraction  (F)  is 
associated  with  the  amount  of  lateral  support  the  graphite 
fibers  receive  from  the  combustion  products.  Using  this 
assumption  the  following  equation  was  used  for  fiber  buckling 
criteria, 
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Ofcr  “ 


F 


(3.5) 


Vf(l-  Vf) 


where  cr' f^r  reduced  buckling  critical  stress  and  F  is 
the  combustion  product  fraction,  that  is,  the  fraction  of  the 
original  stress  due  to  combustion  product  lateral  support. 
Failure  of  the  graphite  fibers  due  to  exceeding  the  reduced 
buckling  stress  as  calculated  above  is  used  by  the  stress 
analysis  program.  Different  values  of  the  combustion  product 
fraction  (F)  are  assumed  to  assess  the  effect  of  this  fraction 
on  failure. 

E.  FIBER  CONSUMPTION 

Fiber  consumption  was  considered  as  the  last  failure  mode. 
As  the  plate  burns,  under  some  conditions,  the  supply  of 
oxygen  in  the  interior  of  the  plate  becomes  depleted.  At  that 
time  the  burning  process  moves  to  the  X=0  surface  where  the 
air  is  being  drawn  into  the  plate.  From  then  on  the  plate 
continues  to  burn  at  the  surface  of  the  plate  slowly 
decreasing  in  thickness  as  the  temperature  profile  across  the 
plate  tends  to  a  constant  value.  This  mode  of  burning  is 
surface  recession. 

The  relatively  constant  temperature  profile  produced 
across  the  plate  in  the  surface  recession  mode,  could  reduce 
the  stresses  seen  by  the  fibers.  If  these  stresses  are  small 
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in  magnitude,  the  possibility  of  fiber  failure  by  consumption 
is  possible.  Surface  recession  fiber  consumption  was 
incorporated  into  the  stress  analysis  code  to  eliminate  burnt 
away  fibers  from  the  analysis. 
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IV.  COMPUTER  PROGRAMS 


A.  MODIFICATIONS  TO  THE  COMBUSTION  CODE 

Only  superficial  modifications  were  made  to  Vatikiotis' 
combustion  code  for  this  study.  First,  the  computer  code  was 
brought  up  to  Fortran  77  language  by  changing  a  few  read  and 
write  format  statements.  Second,  a  single  statement  reading 
numerous  input  variables  from  the  input  file  was  broken  up 
into  individual  read  statements  in  order  to  facilitate  an  easy 
to  read  comment  type  input  file  format.  This  input  file 
containing  a  short  description  of  each  input  variable  is 
contained  in  Appendix  D.  Third,  a  few  output  statements  were 
added  to  facilitate  building  output  data  files  for  use  in  the 
stress  analysis.  The  fourth  and  final  change  to  this  program 
was  to  increase  the  number  of  nodal  points  so  that  the  present 
progrcim  could  handle  up  to  101  nodal  points  vice  the  older  31 
nodal  points.  Increases  in  computer  speed  along  with 
significant  memory  resources  over  the  past  few  years 
facilitated  this  change.  The  amount  of  accuracy  was 
significantly  improved  by  using  a  larger  number  of  nodal 
points,  as  seen  in  both  Figure  8  and  Table  1.  Both  display 
the  same  typical  run  of  the  combustion  code  varying  only  the 
number  of  nodal  points  used.  Temperatures  at  the  locations  of 
X*0"  (the  fire  surface),  X»0.25"  (half  the  medium  thickness). 
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and  X-0.5"  (the  surface  opposite  the  fire)  are  con^ared.  As 
seen,  previous  models  of  the  program  were  not  grid  independent 
for  the  X»0"  location. 
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TABLE  1:  INCREMENTAL  TEMPERATURE  CHANGE  VERSUS  NODAL  POINTS 


USED. 


TIME  OF  RUN:  60  seconds 

HEAT  FLUX  APPLIED:  10,000  BTU/ft^  hr 

DIFFERENTIAL  PRESSURE  ACROSS  TEE  MEDIUM:  0.8  psi 

MEDIUM  THICKNESS:  0.5  inches 

NUMBER 

OF 

NODAL 

POINTS 

USED 

GRAPHITE  TEMPERATURE  °F  AT  LOCATION 

XsQ  inch 

XaO.25  Inch 

X-0.5  Inch 

30 

822.46 

607.35 

486.48 

31 

860.57 

625.83 

486.23 

INCREMENTAL 

CHANGE 

4.4% 

3.0% 

0.05% 

70 

1066.23 

662.77 

493.31 

71 

1068.11 

662.88 

493.14 

INCREMENTAL 

CHANGE 

0.2% 

0.02% 

0.03% 

100 

1099.51 

672.21 

494.42 

101 

1100.33 

672.43 

494.54 

INCREMENTAL 

CHANGE 

0.07% 

0.03% 

0.02% 
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B.  STRESS  ANALYSIS  PROGRAM 


The  stress  program,  as  shown  in  Figure  8,  obtains 
temperature  information  from  an  output  data  file  produced  by 
Vatikiotis'  combustion  code.  Location  data  from  the 
combustion  code  was  presented  in  normalized  (X/L)  form  since 
the  equations  Vatikiotis  used  were  non-dimensionalized.  In 
the  stress  code  these  locations  were  converted  into  absolute 
location  points  (X)  .  Until  burn  off  of  the  fibers  starts 
(i.e.  prior  to  surface  recession)  the  non  dimensionalized 
locations  (X/L)  remain  fixed.  However,  when  surface  recession 
starts,  combustion  moves  to  the  X=0  surface  and  the  thickness 
of  the  material  decreases  with  time.  As  the  thickness  of  the 
material  changes  the  stress  code  continually  verifies  stress 
code  fiber  locations  removing  burned  away  fibers  from  the 
stress  analysis. 

Vatikiotis  in  his  combustion  code  formulation  used  a 
method  of  "smearing"  to  account  for  both  the  fiber  and  air 
filled  void  as  shown  in  Figure  9.  Denoting  the  porosity  of 
the  medium  as  p,  we  have  dV3=pdV  and  dVg={l-p)dV  where  dV  is 
the  total  differential  volume  of  air  and  graphite  fiber,  dV^^ 
is  the  differential  volume  of  air  and  dVg  is  the  differential 
volume  of  graphite.  The  effect  of  this  smearing  process  is  to 
replace  individual  fibers  and  voids,  which  are  extremely  large 
in  number,  by  a  homogenous  material.  This  process  is  utilized 
in  the  basic  differential  equations  presented  in  Appendix  A. 
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STRESS  ANALYSIS  CODE  FLOW  CHART 


SET;  NUMBER  OF  ELEMENTS 
NUMBER  OF  NODES 
L  MATERIAL  PROPERTIES 


COMBUSTION  CODE 
OUTPUT  FILE 


IMSL 

SUBROUTINE 

DLSARO 


SUBROUTINE 
.  READER 


SUBROUTINE 

INTERP 


SUROUTINE  MATRIX^ 
CALCULATE  ELEMENTAL 
^MATRIX  AND  VECTOR  ^ 


READ  IN  COMBUSTION 
CODE  PROFILES 


SUBROUTINE 

MOD2FL 


"  IMSL 
SUBROUTINE 
OLSARO 


INTERPOLATE 
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Figure  8:  Flow  chart  of  the  Stress  Analysis  Code. 
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Figure  9:  The  process  of  "smearing"  used  in  the  combustion 
code. 


The  stress  code  accommodated  the  idea  of  a  homogeneous 
material  and  utilized  temperature  output  by  the  following 
procedure  (see  Figure  10)  .  First  the  graphite  temperatures 
along  with  nodal  point  locations  are  determined  from  the 
combustion  code.  These  combustion  code  locations  are  compared 
to  the  center  location  of  each  ply.  When  the  ply  center 
location  is  bracketed  by  the  combustion  code  locations  the 
stress  code  performs  quadratic  interpolation  to  determine  the 
graphite  fiber  temperature  for  that  ply.  By  this  method, 
accurate  temperature  profiles  can  be  obtained  from  the 
combustion  code  using  101  nodal  points  while  at  the  same  time 
the  nximber  fibers  used  in  the  stress  code  may  be  varied. 
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THREE  POINT  INTERPOLATION  OF  TEMPERATURES 


temperatures  to  stress  code  temperatures. 
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COMBUSTION  CODE  GRAPHITE  TEMPERATURE  STRESS  CODE 


Elemental  matrix  equations  are  determined  using 
Equation  (2.8).  These  elemental  matrices  are  assembled  into 
a  global  set  of  equations  of  the  general  form: 
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where : 

K  =  the  global  stiffness  matrix 
U  =  the  global  displacement  vector 
F  =  the  global  force  vector 
m  =  the  total  number  of  nodal  points  (2n) 
j  and  k  =  are  arbitrary  integers  where  l<j  and  k<m. 
The  boundary  conditions  imposed  on  this  matrix  equation  (4.1)  , 
fit  two  different  forms;  Type  I  and  Type  II.  Their 
application  is  performed  as  follows.  The  Type  I  boundary 
conditions  consist  of  those  that  fit  the  form: 

Uj  =  b  (4.2) 

where : 

Uj  =  displacement  at  node  j , 
b  =  a  known  constant - 
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This  type  of  boundary  condition  is  directly  comparable  to  the 
fixed  end  of  the  fibers  (i.e.  the  locations  of  all  the  odd 
numbered  nodal  points)  where  b=0.  Applying  a  single  boundary 
condition,  Uj=b  to  the  equation  set  (4.1)  above  gives  the 
following  modified  matrix  equation: 
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The  stress  code  was  written  to  apply  each  fixed  end  boundary 
condition  using  the  same  method  as  (4.3)  above. 

The  Type  II  boundary  conditions  fit  the  form  of: 


=  c 


where : 

Uj  and  Uj^  *  are  two  degrees  of  freedom 
a  =  known  constant 
c  =  known  constant. 


(4.4) 
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Application  of  the  general  boundary  condition  of  Equation (4 . 4) 
to  the  matrix  Equation  (4.1)  above  yields: 

(4.5) 
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The  detail  for  developing  this  modified  matrix  equation  is 
given  by  Akin  [Ref.  8] . 

The  Type  II  boundary  condition  of  equation  (4.2)  above  is 
directly  relatable  to  the  fiber  ends  that  move  in  unison, 

C/j  +  (-1)  U,  =  0 

+  (-1)  LTg  =  0 

dg  +  (-1)  Ug  =  0  (4.6) 

U2n-2*  (-1)  =  0 

with  a=-l  and  c=0.  The  Type  II  boundary  conditions  of 
Equation  (4.6)  were  applied  in  the  stress  analysis  program. 

The  final  system  of  equations  with  both  sets  of  boundary 
conditions  applied  was  solved  for  the  displacement  vector  U 
using  the  IMSL  library  subroutine  DLSARG. 

Once  nodal  point  displacements  Uj ( j=l, 2 . . . , 2n)  were  known 
fiber  stresses  cr  were  determined  using  equation  (2.14).  This 
process  was  continued  for  each  time  step  of  output  data 
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obtained  from  the  combustion  code.  The  fiber  stress  analysis 
code  is  included  in  Appendix  B. 
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V.  CASE  STUDIES 


The  combustion  code  requires  input  data  to  start  the 
analysis.  Several  input  data  files  were  produced  encompassing 
various  combinations  of  composite  thickness,  heat  flux, 
differential  pressure  and  combustion  product  fraction  (F)  . 
Specific  values  of  these  variable  parameters  as  used  in  the 
case  studies  are  listed  in  Table  2.  Several  other  important 
parameters  were  treated  as  constants  for  these  studies  as 
listed  in  Table  3.  The  input  data  files  required  by  the 
combustion  code  were  built  covering  all  1008  combinations  of 
the  variable  parameters.  These  data  files  were  processed  by 
the  combustion  code  and  then  the  temperature  profile  data  was 
further  processed  for  stress  failure  information.  Most  of  the 
lower  heat  flux  results,  primarily  below  4000  btu/  hr-ft^, 
produced  relatively  minor  increases  in  plate  temperatures 
which  resulted  in  insignificant  stresses  and  no  failure.  But 
the  majority  of  runs  gave  interesting  results.  The  final 
results  of  this  effort  are  provided  in  graphical  form  in  the 
following  case  studies. 
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A.  TEMPERATURE  PROFILES 


As  Stated  in  previously,  induced  stresses  are  dependent 
upon  the  fiber  temperatures.  Figures  11  through  16  show  the 
temperature  profiles  of  some  representative  combustion  code 
runs  up  to  the  point  of  initial  fiber  failure.  Note  that  the 
maximum  temperature  is  always  at  the  bottom  of  the  plate  where 
the  heat  flux  is  being  applied  and  the  air  is  coming  in,  while 
the  cooler  fibers  are  at  the  opposite  surface.  The  main 
concept  to  bare  in  mind  for  these  studies  is  that  graphite 
fibers  contract  upon  heating,  so  the  hotter  fibers  will  be  in 
a  tensile  mode  (see  Figure  4) .  Since  the  tensile  strength  of 
these  hot  fibers  is  very  high,  these  fibers  will  not  fail. 
However,  the  fibers  that  are  cooler  will  be  in  a  compression 
mode  with  a  low  critical  buckling  strength  and  they  will  fail 
by  buckling. 

A  comparison  of  Figures  11,  12,  and  13  show  the  effect  of 
increasing  the  heat  flux  on  a  constant  thickness  material. 
Notice  that  for  the  lower  heat  flvixes,  the  temperatures  are 
lower  and  the  profile  is  more  gradual.  This,  of  course,  is 
due  to  the  lower  heat  flux  applied  but  also  due  to  the  ability 
of  the  air  flow  to  convectively  transfer  the  heat  away  from 
the  heated  surface.  If  the  heat  flux  is  low,  the  flow  of  air 
through  the  medium  has  more  time  to  distribute  the  heat 
through  the  material  before  failure  occurs.  If  the  heat  flux 
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is  high,  the  surface  temperature  rises  quickly  and  heating  is 
confined  primarily  to  the  area  close  to  the  heated  surface. 

A  comparison  of  Figures  11  and  16  indicates  the  changes 
in  the  temperature  profile  when  the  highest  and  lowest 
differential  pressures  are  applied  for  a  given  heat  flux  and 
material  thickness.  As  stated  previously,  the  convective  heat 
transfer  in  the  medium  has  a  large  effect  on  the  temperature 
profiles  for  the  lower  heat  fluxes.  For  a  constant  thickness, 
the  pressure  gradient  AP/L  will  be  higher  for  a  higher 
differential  pressure.  This  higher  pressure  gradient  causes 
a  higher  air  flow  through  the  material  which  increases  the 
convective  heat  transfer.  This  increased  heat  transfer  causes 
both  a  lower  and  more  gradual  temperature  profile  in  the  low 
heat  flux  cases.  However,  in  the  high  heat  flux  cases  this 
effect  is  insignificant  as  seen  in  Figures  12  and  14  for  a 
heat  flux  of  10,000  btu/hr-ft^. 

As  a  final  comparison  Figures  12  and  15  are  compared. 
These  figures  show  the  dependence  of  the  temperature  profile 
on  thickness  of  the  material  while  holding  differential 
pressure  constant.  The  temperatures  at  the  surface  of  the 
1.0"  plate  are  higher  than  the  0.5"  plate.  This  can  be 
explained  as  follows.  If  the  AP  is  a  fixed  value,  then  the 
thicker  plate  has  a  smaller  pressure  gradient  (AP/L)  because 
the  thickness  L  is  larger.  This  smaller  AP/L  gradient  will 
cause  less  air  flow  through  the  thicker  material  resulting  in 
a  lower  cunount  of  convective  heat  transfer.  The  end  result  for 
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the  thicker  material  is  a  higher  and  more  steep  temperature 
profile.  This  profile  causes  more  fibers  to  be  in  the 
compressive  mode  and  subsequently  fail  by  buckling. 
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TEMPERATURE  PROFILE 

WITH  RESPECT  TO  TIME 


TIME  TO  RRST  FAILURE-109.2  sec 
HEAT  FLUX-3,000  BTU/HR  FT* 
0.5  INCH  THICKNESS 
I  AP-0.7PSI 


Slgur*  11:  Temperacure  profile  up  to  first  failure  for  0.5", 
3,000  btu/br  ft^  and  0.7  psi. 
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Figur*  12:  Ten5)erature  profile  up  to  first  failure  for  0.5", 
10,000  btu/hr  ft^  and  0.7  psi. 
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TEMPERATURE  PRORLE 

WITH  RESPECT  TO  TIME 

TIME  TO  INITAL  FAILURE=4.9  sec 
HEAT  R.UX=25,000  BTU/HR  FT^ 
0.5  INCH  THICKNESS 


Figure  13:  Temperature  profile  up  to  first  failure  for  0.5", 
25,000  btu/hr  ft^  and  0.7  psi. 
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TEMPERATURE  PROFILE 

WITH  RESPECT  TO  TIME 


TIME  TO  FIRST  FAILURE-12.0  sec 
HEAT  FLUX-10,000  BTU/HR  FT^ 
0.5  INCH  THICKNESS 
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Flgur*  14:  Temperature  profile  up  to  first  failure  for  0.5 
10,000  btu/hr  ft^  and  0.05  psi. 
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TEMPERATURE  PROFILE 

WITH  RESPECT  TO  TIME 


TIME  TO  FIRST  FAILURE-27.3  sec 
HEAT  FLUX-10,000  BTU/HR  FT* 
1.0  INCH  THICKNESS 


Flgur*  15:  Temperature  profile  up  to  first  failure  for  1.0", 
10,000  btu/hr  ft^  and  0.7  psi. 
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TEMPERATURE  PROFILE 

WITH  RESPECT  TO  TIME 


TIME  TO  FIRST  FAILURE=58.3  sec 
HEAT  FLUX=3,000  BTU/HR  FT^ 
0.5  INCH  THICKNESS 
AP=0.05  PSI 


Figure  16;  Temperature  profile  up  to  first  failure  for  0.5", 
3,000  btu/hr  ft^  and  0.05  psi. 


47 


B.  TIME  TO  INITIAL  FAILURE  (AP  EFFECT) 

Figures  17  through  20  present  the  time  to  initial  failure 
showing  its  dependency  on  both  heat  flux  and  differential 
pressure.  As  heat  flux  decreases,  the  time  to  initial  failure 
progressively  increases.  For  each  specific  AP  there  appears 
to  be  a  lower  bound  of  heat  flux  below  which  failure  will 
never  occur.  The  figures  also  show  that  this  lower  bound 
value  increases  with  increasing  AP.  However,  this  lower  bound 
cannot  be  exactly  determined  through  use  of  a  computer  code. 
Through  the  use  of  finite  increments  of  heat  flux,  given 
previously  in  Table  2,  this  lower  bound  for  failure  was 
roughly  determined.  As  seen  in  Figures  17  through  20  this 
lower  bound  of  heat  flux  was  found  to  within  500  btu/hr-ft^. 

For  all  thicknesses,  the  effect  of  pressure  is  relatively 
small  except  at  the  lower  heat  fluxes  where  the  effect  of  air 
flow  through  the  material  has  a  large  effect  on  the  convective 
heat  transfer.  These  pressure  effects  become  ineffectual  at 
higher  heat  fluxes. 
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TIME  TO  INITIAL  FIBER  FAILURE 
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Figure  18:  The  time  to  initial  failure  (effect  of 
differential  pressure)  for  0.75",  and  F«0.5%. 
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Figure  20:  The  time  to  initial  failure  (effect  of 
differential  pressure)  for  0.25",  and  F»0.5%. 
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C.  TIME  TO  INITIAL  FAILURE  (P  EFFECT) 

Figures  21  through  24  demonstrate  how  the  critical 
buckling  strength,  and  hence  time  to  initial  failure,  varies 
with  the  combustion  product  factor  (F)  .  The  value  of  F 
greatly  effects  not  only  the  time  to  failure  but  determines 
whether  or  not  failure  will  occur.  As  F  increases  there  is 
more  lateral  support  for  the  fibers  so  the  critical  buckling 
strength  increases.  This  increase  results  in  fewer  fibers 
failing  by  buckling.  In  the  case  of  the  0.25"  plate  none  of 
the  fibers  reached  the  critical  buckling  strength  when  F=1.0%. 
All  figures  show  that  for  a  given  heat  flux,  as  F  increases, 
the  time  to  initial  failure  increases. 

The  result  of  this  analysis  shows  that  the  value  assumed 
for  the  combustion  product  fraction  (F)  greatly  influences  the 
amount  of  time  to  initial  failure  of  the  material.  The 
figures  show  that  when  F  increases  the  resistance  to  buckling 
increases  and  buckling  is  less  likely  to  occur. 
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Flgur«  21:  The  time  to  initial  failure  (effect  of  combustion 
product  fraction)  for  1.0",  and  AP-0.7  psi. 


Figure  22:  The  time  to  initial  failure  (effect  of  combustion 


product  fraction)  for  0.75",  and  AP-0.7  psi. 
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TIME  TO  INITIAL  FIBER  FAILURE 

(EFFECT  OF  F) 
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product  fraction)  for  0.5",  and  AP-0.7  psi. 
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Figure  24:  The  time  to  initial  failure  (effect  of  combustion 
product  fraction)  for  0.25",  and  AP-0.7  psi. 
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O.  PROGRESSION  OF  MATERIAL  FAILURE  (AP  EFFECT) 

Figures  25  and  26  track  the  percentage  of  failed  fibers 
from  the  time  of  the  first  observed  fiber  failures.  Due  to 
the  0.5  second  time  step  used  for  temperature  profile  output, 
some  of  these  curves  start  at  time  zero  with  a  significant 
number  failed  fibers.  Actually  the  failure  represented  by 
these  curves  commenced  some  time  between  the  time  of  zero  of 
the  figures  and  the  previous  0.5  second  interval.  Since 
observation  of  trends  in  this  case  are  important,  the 
increased  computer  time  to  determine  a  more  precir-'  starting 
time  was  not  justified. 

The  general  trend  is,  the  higher  the  differential 
pressure,  the  lower  the  percentage  of  failure.  This  is 
apparently  is  due  to  the  more  gradual  temperature  profiles 
caused  by  the  higher  air  flows  through  the  material.  For  the 
higher  differential  pressures/higher  air  flowrates,  more 
fibers  are  heated  to  approximately  the  same  temperature,  and 
are  therefore  in  the  tensile  mode  prior  to  the  time  of  initial 
fiber.  These  tensile  fibers  do  not  fail. 

Another  item  worth  mentioning  is  the  rapidity  of  failure. 
The  materials  progressed  to  their  maximum  percentage  of  fiber 
failure  very  rapidly.  For  the  1.0"  thicker  material  this 
maximum  failure  percentage  is  reached  within  3  seconds  of  the 
first  fiber  failure  and  for  the  0.5"  thinner  material 
approximately  4  seconds.  These  times  are  independent  of 
differential  pressures. 
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TIME  FROM  FIRST  FAILURE  (sec) 


Figure  25:  The  progression  of  material  failure  (the  effect 
of  AP)  for  0.5",  and  5,000  btu/  hr-ft^. 


TIME  FROM  RRST  FAILURE  (sec) 


Figure  26:  The  progression  of  material  failure  (the  effect 
of  AP)  for  1.0",  and  5,000  btu/hr-ft^. 


55 


E.  PROGRESSION  OF  MATERIAL  FAILURE  (EFFECT  OF  HEAT  FLUX) 

Figures  27  and  28  indicate  that  as  heat  flux  increases, 
time  to  maximvim  failure  changes  only  marginally,  effecting  the 
3  to  4  second  failure  time  by  only  1  second. 
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0  1  2  3  4  5 


TIME  FROM  FIRST  FAILURE  (sec) 

Figure  27:  The  progression  of  material  failure  (the  effect 
of  heat  flux)  for  0.5",  and  AP=0.7  psi. 


TIME  FROM  RRST  FAILURE  (sec) 


Figure  28:  The  progression  of  material  failure  (the  effect 
of  heat  flux)  for  1.0",  and  AP-0.7. 
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F.  MAXIMOM  FAILURE  (EFFECT  OF  AP) 

Figures  29  through  32  are  a  consolidation  the  maximum 
fiber  failure  for  material  failure.  These  maximum  percentages 
varied  slightly  with  pressure  at  lower  heat  fluxes.  As  shown, 
for  a  given  material  thickness,  higher  heat  fluxes  give  higher 
percentage  of  failed  fibers.  This  is  due  to  the  more  non- 
uniform  temperature  profile  discussed  previously.  The  maximum 
percentage  failure  seen  in  these  studies  is  the  case  of  the 
one  inch  plate  at  25,000  btu/hr-ft^,  which  reached  a  maximum 
of  97.5%. 
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Figure  29:  The  maximxiin  percentage  of  failure  (the  effect  of 


AP)  for  1.0",  and  F*0.5%. 


Figure  30:  The  maucimum  percentage  of  failure  (the  effect  of 
AP)  for  0.75",  and  F-0.5%. 
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MAXIMUM  FAILURE  VS  HEAT  FLUX 

(EFFECT  OF  AP) 
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Figure  31:  The  maximum  percentage  of  failure  {the  effect  of 
AP)  for  0.5",  and  F=0.5%. 
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Figure  32:  The  maximum  percentage  of  failure  (the  effect  of 
AP)  for  0.25",  and  F-0.5%. 
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G.  MAXIMUM  FAILX7RE  (EFFECT  OF  F) 

The  effect  of  the  combustion  product  fraction  is  shown  in 
Figures  33  and  34.  As  expected,  assuming  a  higher  value  of 
the  combustion  product  fraction  (F)  reduces  the  amount  of 
failure  because  of  the  increase  in  buckling  strength.  For  the 
0.25"  plate  there  was  no  failure  above  a  factor  of  1.0%  for 
range  of  heat  fluxes  investigated. 

These  trends  were  independent  of  differential  pressure. 
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Figure  33;  The  maximum  percentage  of  failure  (the  effect  of 


combustion  product  fraction)  0.5",  and  AP=0.7  psi. 


Figure  34:  The  maximum  percentage  of  failure  (the  effect  of 
combustion  product  fraction)  0.25"  and  AP-0.7. 
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H.  MAXnCDM  FAILURE  (EFFECT  OF  THICKNESS) 

As  seen  in  Figures  35  and  36,  the  thickness  of  the 
material  has  minimal  effect  on  the  final  percentage  of 
material  failure  for  higher  heat  fluxes.  However  at  lower 
heat  fluxes  the  thinner  materials  experience  significantly 
less  failure.  In  all  cases  the  final  percentage  of  fiber 
failure  was  between  42%  and  97.5%. 
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Figure  36:  The  maximum  percentage  of  failure  (the  effect  of 
material  thickness)  F-0.5%,  and  AP~0.05. 
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VI.  CONCLUSIONS 


The  conclusions  drawn  from  this  study  are  as  follows: 

#  The  amount  of  time  required  before  onset  of  fiber  failure 
ranged  from  a  very  small  3  seconds  (for  a  0.25"  plate  at 
25,000  btu/hr-ft^)  up  to  a  maximum  of  210  seconds  (for  a 
1.0"  plate  at  2,000  btu/hr-ft^) .  Given  the  reaction  times 
of  highly  trained  fire  fighting  teams  onboard  aircraft 
carriers,  it  is  quite  conceivable  that  failure  may  be 
averted  for  the  longer  onset  of  failure  times. 

#  Once  the  failure  of  carbon  graphite -epoxy  composite 
material  begins  maximxim  failure  occurs  very  quickly  in 
less  than  4  seconds.  This  time  is  of  such  short  duration 
that  reducing  the  amount  of  maximum  failure  once  failure 
begins  is  not  possible. 

#  The  predominant  factors  effecting  survivability  of  carbon- 
graphite  epoxy  composite  materials  are  the  heat  flux 
applied,  the  combustion  product  factor  and  the  thickness 
of  the  material.  The  differential  pressure  effect  varies 
from  very  important  at  low  heat  fluxes  to  insignificant  at 
high  heat  fluxes.  However,  the  effect  of  differential 
pressure  is  mitigated  by  the  fact  that  response  times  are 
in  all  cases  very  small  at  the  high  heat  fluxes.  Once 
failure  occurs,  maximum  failures  in  the  range  of  80-90% 
preclude  any  future  use  of  the  material . 
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APPENDIX  A:  COMBUSTION  CODE  EQUATIONS 

LIST  OF  EQUATION  SYMBOLS 

H  -  dynamic  viscosity  of  air 
4)  -  oxygen  concentration 

-  air  density 
Pg  -  graphite  density 
Cg^  -  specific  heat  capacity  of  air 
Cg  -  specific  heat  capacity  of  graphite 
Dg  -  Mass  diffusion  coefficient 
h  -  convection  heat  transfer  coefficient 
kg  -  air  thermal  conductivity 

kg  -  effective  thermal  conductivity  of  porus  graphite 

kj.  -  pseudo  thermal  conductivity  to  account  for  radiation 

P  -  pressure 

p  -  porosity 

^02  ■  oxygen  reaction  rate 

Rg  -  graphite  reaction  rate 

Tg  -  air  temperature  (°R) 

Tg  -  graphite  temperature  (°R) 
t  -  time 

u  -  pore  velocity  (of  air) 

z  -  specific  internal  area 
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PINAL  FIELD  EQUATIONS 

Following  are  the  equations  used  in  the  Vatikiotis 
Combustion  code.  These  equations  are  presented  here  in  their 
form  prior  to  non-dimensionalization.  These  equations  were 
non- dimens ionali zed  prior  to  use  in  the  actual  code. 

THE  GRAPHITE  HEAT  TRANSFER  EQUATION 


THE  INTERNAL  AIR  HEAT  TRANSFER  EQUATION 


*  hziT  -T,)  ■> 


=  p 

OX 


dT, 

at 


OXYGEN  TRANSPORT  EQUATION 


44)  -  -i-(pu<i>)  -  R, 

Ox  Ox 


THE  INTERNAL  FLOW  EQUATION 


P  ^  ^  1  ^Pa  ^  ^  -  H  ^(PPa^  ^ 

dx^  Pa  m  dx  dx  dx  p^m  dt 
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BOUNDARY  CONDITIONS  IMPOSED 

For  the  X=0  location, 

=  T. 

<{,  = 

p  ^  p. 


p  =  p^  +  ^p 


where : 

qgt  -  the  heat  flux  applied 
T„  =  ambient  temperature  (°R) 

<t>„  =  ambient  oxygen  concentration 

hj^  and  h2  *  convection  heat  transfer  coefficients 
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INITIAL  CONDITIONS 


There  are  a  number  of  ways  to  start  these  "combustion” 
problems.  One  way  is  to  set  the  oxygen  concentrations  and  the 
air  and  graphite  temperatures  to  ambient  conditions  then  apply 
an  external  heat  source.  This  "natural"  starting  scheme 
requires  some  amount  of  computer  effort  to  bring  the  plate  up 
to  conditions  where  temperatures  have  increased  significantly. 
Instead  of  starting  the  problem  at  ambient  conditions, 
computational  effort  can  be  reduced  by  starting  the  problem 
with  the  plate  at  elevated  temperatures. 

In  this  particular  study  the  initial  conditions  were; 

Graphite  temperature  400°F 

Air  temperature  400°F 

Oxygen  concentration  0.0172  Ib/ft^  (ambient  level). 
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APPENDIX  B:  STRESS  ANALYSIS  CODE 


♦EDWARD  FAXLANGER  JULY  1992 

*  THIS  PROGRAM  CALCULATES  STRESSES  INDUCED  IN  LONGITUDINAL  FIBERS 

*  RESULTING  FROM  TEMPERATURE  VARIATION  ACROSS  THE  MEDIUM 

4rW**ilir1lr4r**'*****ii^il»*irir^1ir4nAr**4r*^1lp****4r*4rik**i*****i^ilr4r4r***#'Ar***i^**4r**1lr*'*r** 

C  DESCRIPTION  OF  VARIABLES 
C  U  =  DISPLACEMENTS  OF  NODAL  POINTS 

C  BC  =  VALUE  OF  THE  TYPE  I  BOUNDARY  CONDITION  FOR  NODAL  POINT  IBC 

C  A  =  MATRIX  CONTAINING  INDIVIDUAL  STIFFNESS 

C  B  =  MATRIX  CONTAINING  INDIVIDUAL  FORCE  VECTOR 

C  AA  =  MATRIX  CONTAINING  GLOBAL  STIFFNESS 

C  BB  =  MATRIX  CONTAINING  GLOBAL  FORCE  VECTOR 

C  Y  =  LOCATIONS  OF  THE  INDIVIDUAL  FIBER  (FROM  BOTTOM  OF  PLATE) 

C  TEMP  =  INTERPOLATED  TEMPERATURE  OF  THE  FIBER 
C  DIAM  =  INTERPOLATED  DIAMETER  OF  THE  FIBER 
C  SIGMA  =  CALCULATED  STRESS  OF  THE  INDIVIDUAL  FIBER 
C  YIN  =  LOCATIONS  OF  NODAL  POINTS  FROM  THE  COMBUSTION  CODE 
C  TEMPIN=TEMPERATURES  FROM  THE  COMBUSTION  CODE 

C  DIAMIN=THE  DIAMETER  AT  THE  NODAL  POINT  FROM  THE  COMBUSTION  CODE 
C  IBC  =  THE  NUMBER  OF  THE  NODE  HAVING  A  TYPE  I  BOUNDARY  CONDITION 

C  NODES®  A  MATRIX  USED  FOR  ELEMENT  CONNECTIVITY  TO  THEIR  NODAL  POINTS 
C  Y1  =  LOCATION  OF  THE  NODAL  POINTS  BROUGHT  IN  FROM  VATIKIOTIS  CODE 
C  NFLAG  »  A  FLAG  TO  INDICATE  TYPE  OF  FAILURE 

C  NELEM  =A  MATRIX  CONTAINING  THE  NAMES  OF  THE  UNBROKEN  FIBERS 
C  NE  ®  NUMBER  OF  FIBER  ELEMENTS 
C  NNODE  =  TOTAL  NUMBER  OF  FIBER  NODES 

C  PERCNT®  PERCENTAGE  OF  ORIGINAL  BUCKLING  STRENGTH  USED 
C  TSTREN®  TENSILE  STRENGTH  OF  FIBERS 

C  CSTREN®  COMPRESSIVE  STRENGTH  OF  FIBERS  (NO  BUCKLING) 

C  EFIB  =  ELASTIC  MODULUS  OF  THE  FIBERS 
C  VF  =  FIBER  VOLUME  FRACTION  (ASSUMED  CONSTANT) 

C  EMAT  =  ELASTIC  ITODULUS  OF  THE  MATRIX 
C  POISM  =  POISON  RATIO  OF  THE  MATRIX 
C  ALPHA  ®  THERMAL  EXPANSION  COEFFICIENT  (AXIAL) 

C  WIDTH  ®  THE  ASSUMED  LENGTH  OF  EACH  FIBER  (ASSUMED  TO  BE  ONE) 

C  TEMREF®  REFERENCE  TEMPERATURE  (WHERE  STRAIN  EQUAL  ZERO) 

C  NNP  =  THE  NUMBER  OF  NODAL  POINTS  READ  IN  FROM  COMBUSTION  CODE 

C  TIME  ®  TIME  OF  THE  TEMPERATURE  PROFILE  BEING  ANALYZED 
C  THICK  =  THE  ORIGINAL  THICKNESS  OF  THE  MATERIAL  (BEFORE  BURNING) 

C  HEIGHT®  THE  THICKNESS  OF  THE  MATERIAL  AT  THE  PRESENT  TIME 
C  NEUSED®  THE  NUMBER  OF  UNBROKEN  FIBERS  AT  PRESENT  TIME 
C  NELEM  ®  A  STORAGE  OF  UNBROKEN  FIBERS  NAMES 
C  GMAT  ®  SHEAR  MODULUS  OF  THE  EPOXY  MATRIX  MATERIAL 
C 

C  MAIN  PROGRAM 

IMPLICIT  REAL*8 (A-H,0-Z) 

C  THE  FOLLOWING  PARAMETER  N1  SETS  THE  DIMENSION  OF  ARRAYS  AND 
C  NEEDS  TO  BE  EQUAL  TO  THE  HIGHER  OF  1)THE  NUMBER  OF  NODES  BEING 
C  READ  IN  ♦♦OR**  2)  THE  NUMBER  OF  ELEMENTS  BEING  ANALYZED 
PARAMETER  (Nl®101) 

PARAMETER  (NN®2*N1) 

C  COMMON  STATEMENT  FOR  MATERIAL  PROPERTIES 
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COMMON/PROPER/TSTREN, CSTREN, EFIB , EMAT, ALPHA, WIDTH, THICK, POISM, VF 
C 

REAL*8  U(NN) ,BC(NN) ,A(2,2) ,B(2) ,BB(NN) ,AA(NN,NN) ,Y(N1) , 
&TEMP(N1) ,DIAM(N1) ,SIGMA(N1) ,YIN(N1) ,TEMPIN(N1) ,DIAMIN(N1) 
DIMENSION  IBC(NN) .NODES (NN,NN) , INDEX(2) ,NFLAG(N1) ,NELEM(N1) 

♦THE  FOLLOWING  3  STATEMENTS  INCREASE  THE  NUMBER  OF  EQUATIONS 
♦THAT  THE  EQUATION  SOLVER  DLSARG  WILL  HANDLE.  IF  LARGER  SPACE  IS 

♦  NEEDED  DLSARG  WILL  GIVE  AN  ERROR  MESSAGE  AND  RECOMMEND  LAR3ER 

♦  WORK  SPACES 

COMMON/WORKSP/RWKSP 
REAL  RWKSP (82236) 

CALL  IWKIN(82236) 

C  DATA  FILE  FOR  OUTPUT  OF  DATA 

OPEN(UNIT=ll,  FILE=' /TEST31  STRESS  Bl' , STATUS='OLD' ) 

C  DATA  FILE  FOR  READING  INPUT  DATA 

OPEN(UNIT=12,  FILE=' /TEST31  TEMPER  Bl' , STATUS =' OLD' ) 

C  NE  IS  THE  NUMBER  OF  FIBERS  ASSUMED  FOR  ANALYSIS 
NE=30 

NNODE=NE^2 

C  ENTER  THE  PERCENTAGE  OF  MATRIX  SUPPORT 
PERCNT=.005 

C  FIBER  CHARACTERISTICS;  LONGITUDINAL  COMPRESSIVE  AND  TENSILE 
C  STRENGTHS  (PSD ,  ELASTIC  MODULUS  AND  THERMAL  COEFFICIENT 
TSTREN=300.0D3 
CSTREN=-260.0D3 
EFIB=31.0D6 
EMAT= . 5D6 
VF=.196 
POISM=.35 
ALPHA»-0.5SD-6 

C  WIDTH  OF  SIMULATED  HOLE  (LENGTH  OF  FIBERS)  ASSUMED  TO  BE  UNITY 
C  SINCE  IT  HAS  NO  EFFECT  ON  PROBLEM 
WIDTHsl . ODO 

C  REFERENCE  TEMPERATURE  IN  DEGREES  F  WHERE  INITIAL  STRESSES  WERE  ZERO 
TEMREF=  60. DO 

C  INITIALIZE  FLAGS  TO  ZERO  (NFLAG  KEEPS  TRACK  OF  BROKEN  FIBERS) 

C  (NSTOP  STOPS  THE  PROGRAM  ) 

NSTOP=0 
DO  M=1,NE 
NFLAG (M)  =0 
END  DO 


♦  READ  IN  OF  INFORMATION  FROM  COMBUSTION  CODE  OUTPUT 
READ(12,«)NNP 
CALL  CHECK (NNP, NSTOP) 

IF (NSTOP. NE.O) GOTO  100 

CALL  READER (N1 , NNP , TIME , HEIGHT, YIN, TEMPIN, DIAMIN, NSTOP) 
IF (NSTOP. NE.O) GOTO  100 

♦LOOP  FOR  FILLING  ALL  ELEMENT  HEIGHTS  ASSUMING  EVEN  SPACING 
THICK=HEIGHT 
YINC=THICK/ (NE) 

YSTART=YINC/2 
DO  Msl,NE 

Y (M) =YSTART+YINC* (M- 1 ) 

END  DO 


71 


o  n  n  n  o 


200  BURN=THICK-HEIGHT 
C  INITIALIZE  ARRAYS  TO  ZERO 
DO  1=1, NN 
BB(I) =0.0D0 
DO  J=1,NN 
AA(I, J) =0.0D0 
END  DO 
END  DO 

DO  1=1, NE 

C  IF  THE  FIBERS  ARE  ALREADY  CONSUMED  BY  COMBUSTION  (OUTSIDE  EDGE) 

C  NULL  THEIR  EFFECT  BY  SETTING  THEIR  FLAG  TO  1 
IF(Y(I) .LT. BURN) THEN 
NFLAG(I) =1 
END  IF 

DETERMINE  THE  TEMPERATURE  AND  DIAMETER  OF  NON- FAILED  FIBERS  THIS 
LOOP  IS  FOR  INTERPOLATION  OF  DATA  RECEIVED  FROM  "STRESS  INPUT"  FILE. 
THE  RECEIVED  DATA  IS  IN  RELATIVE  SPACING  AND  MUST  BE  CHANGED  TO 
ABSOLUTE  SPACING  IN  ORDER  TO  DETERMINE  PERCENTAGE  FAILURE  OF  FIBERS 
INTERPOLATE  NON-BROKEN  FIBERS  ONLY 
IF(NFLAG(I) .EQ.O)  THEN 
DO  N=l,NNP-2 

IF( (Y(I) .GE.YIN(N) +BURN) -AND. (Y(I) . LE . YIN (N+2) +BURN) ) THEN 
CALL  INTERP (YIN (N) +BORN, YIN (N+1) +BURN, YIN (N+2 ) +BURN, 

&DIAMIN(N)  ,DIAMIN(N+1)  ,DIAMIN(N+2)  ,Y(I)  ,DIAM(I)  ) 

CALL  INTERP (YIN (N) +BURN, YIN (N+1) +BURN, YIN (N+2 ) +BURN, 

&TEMPIN(N) ,TEMPIN(N+1) ,TEMPIN(N+2) ,Y(I) , TEMP (I) ) 

END  IF 
END  DO 
END  IF 
END  DO 

C  LOOP  FOR  DETERMINING  THE  ELEMENTS  THAT  ARE  NOT  BROKEN  (IE  THE 
C  ONES  THAT  ARE  TO  BE  ANALIZED) 

IJEUSED=0 
DO  1=1, NE 

IF(NFLAG(I) .EQ.O)  THEN 
NEUSED=NEUSED+1 
NELEM(NEUSED) =I 
END  IF 
END  DO 

IF (NEUSED . LT . 2 ) THEN 

WRITE (11,*) 'ALL  BUT  TWO  ELEMENTS  HAVE  BROKEN' 

WRITE (11,*) 'ANALYSIS  IS  STOPPED /PROGRAM  STOPPED' 

GOTO  100 
END  IF 

NNUSED-NEUSED*2 

*  LOOP  FOR  ESTABLISHING  CONNECTIVITY  OF  NODES  TO  THEIR  ELEMENTS 

*  FOR  USE  IN  BUILDING  GLOBAL  MATRIX  AND  GLOBAL  VECTOR 

L=1 

DO  1=1, NEUSED 
DO  J=l,2 
NODES ( I, J) =L 
L=L+1 
END  DO 
END  DO 


*  L(X)P  TO  BUILD  GLOBAL  MATRIX (AA)  AND  VECTOR (BB) 
DO  IE=1, NEUSED 
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CALL  MATRIX  (A, B,  TEMP  (NELEMdE)  )  -TEMREF,  DIAM (NELEM ( IE)  )  ) 

INDEX ( 1 ) =NODES ( IE , 1 ) 

INDEX ( 2 ) =NODES ( IE , 2 ) 

DO  1=1,2 
II=INDEX(I) 

BB(II) =BB(II) +B(I) 

DO  J=l,2 
JJ=INDEX(J) 

AA(II, JJ) =AA{II, JJ) +A(I, J) 

END  DO 
END  DO 
END  DO 

*LOOP  FOR  TYPE  I  BOUNDARY  CONDITION  INITIALIZATION 
DO  I=1,NNUSED-1, 2 
IBC(I) =I 
BC(I) =0.0D0 
END  DO 

*  APPLY  TYPE  1  BOUNDARY  CONDITIONS 
DO  I=1,NNUSED-1, 2 
J=IBC(I) 

DO  K=1,NNUSED 

BB(K)  =BB(K)  -BC(I)  *AA(K,  J) 

AA(K, J) =0.0 
AA(J,K) =0.0 
AA(J, J) =1.0 
BB(J)=BC(I) 

END  DO 
END  DO 

C  APPLY  TYPE  2  BOUNDARY  CONDITIONS 
DO  J=2,NNUSED-2,2 

CALL  MOD2FL(NN,NNUSED, J,J+2, -1.0D0,0.0D0,AA,BB) 

END  DO 

C  SOLVE  THE  SYSTEM  OF  EQUATIONS  USING  IMSL  LIBRARY 
CALL  DLSARG(NNUSED,AA,NN,BB,1,U) 

DO  1=1, NE 
SIGMA (I) =0.0D0 
END  DO 

DO  I=1,NEUSED 
STRAIN=U(I*2) /WIDTH 

SIGMA (NELEM (I) ) =EFIB* (STRAIN- ALPHA* (TEMP (NELEM (I) ) -TEMREF) ) 

END  DO 

C  LOOP  TO  CALCULATE  STRESSES  OF  EACH  ELEMENT  AND  OUTPUT  RESULTS 
C  TO  OUTPUT  FILE 

WRITE (11, 9) TIME 
WRITE (11, 6) 

DO  1=1, NE 

IF(NFLAG(I) .NE.0)WRITE(11,8)  1,NFLAG(I) 

IF(NPLAG(I) .EQ.O) WRITE (11, 7) I. SIGMA (I) , DIAM (I) 

END  DO 

6  FORMAT (5X, 'FIBER  NUMBER' ,5X, 'FIBER  STRESS  (PSD ', 2X, ' DIAMETER' ) 

7  FORMAT(12X,I3,11X,G10.4,6X,G10,4,3X,G11.4,2X,I1) 

8  FORMAT (17X, I 3, 4X, 'PREVIOUSLY  BROKEN,  NFLAG=',I1) 

9  FORMATdX, 'TIME=' ,F8.3,  '  SECONDS') 

CALL  BROKEN (N1 , NE , SIGMA, DIAM, NFLAG, BROKE , PERCNT) 

WRITE (11, *) ' *****************************************' 
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oooonn  n n  on  on 


WRITE (11, *) 

CALL  READER (N1 , NNP, TIME, HEIGHT, YIN, TEMPIN, DIAMIN, NSTOP) 
IF(NSTOP.NE.O)GOTO  100 
GOTO  200 
100  END 


SUBROUTINE  MATRIX ( A, B , TEMP , DIAM) 

THIS  SUBROUTINE  CALCULATES  THE  ELEMENT  LEVEL 
STIFFNESS  MATRIX (A)  AND  VECTOR (B) 

IMPLICIT  RKAL*8 (A-H,0-Z) 

COMMON/PROPER/TSTREN, CSTREN, EFIB, EMAT, ALPHA, WIDTH, THICK, POISM, VF 
REAL*8  A{2,2) ,B(2) 

PI=4*ATAN(1 .0) 

AREA=DIAM*DIAM*PI/4 . 0 

S=AREA*EFIB/WIDTH 

A(l,l) =S 

A(l,2) =-S 

A(2,l) =-S 

A(2,2) =S 

B (2) =EFIB*AREA*ALPHA*TEMP 
B(l)  =-B(2) 

END 


SUBROUTINE  BROKEN (N1 , NE , STRESS , DIAM, NFLAG, BROKE , PERCNT) 

THIS  SUBROUTINE  DETERMINES  THE  NUMBER  OF  BROKEN  FIBERS  (BROKE) 

DUE  TO  TENSILE  AND  COMPOSITE  BUCKLING 
NFLAG  EQUALS  0 : FOR  UNFAILED  FIBER 

l:FOR  BURNT  AWAY  FIBERS 
2: FOR  OVER  STRESS  BY  TENSION 
3: FOR  OVER  STRESS  BY  BUCKLING 
IMPLICIT  REAL*8 (A-H,0-Z) 

COMMON/PROPER/TSTREN, CSTREN, EFIB, EMAT, ALPHA, WIDTH, THICK, POISM, VF 
REAL*8  STRESS (Nl) , DIAM (Nl) 

DIMENSION  NFLAG (Nl) 

PI=4*ATAN(1.0) 

C  CALCULATE  COMPOSITE  BUCKLING  STRENGTH  OF  THE  FIBER 
C  AS  A  FRACTION  OF  THE  ORIGINAL  COMPOSITE  BUCKLING  STRENGTH 
GMAT=EMAT/ (2 . * (1 . +POISM)  ) 

BUCKLE=-GMAT/ (1 . -VF) *PERCNT 
DO  I:rl,NE 

C  DO  NOT  CHECK  PREVIOUSLY  BROKEN  FIBERS 
IF(NFLAG(I) .EQ.O)  THEN 
C  CHECK  FOR  TENSILE  STRESS  FAILURE 

IF(STRESS (I) .GE.TSTREN)  THEN 
NFLAG (I) =2 
END  IF 

C  CHECK  FOR  BUCKLING  FAILURE 

IF (STRESS (I) .LE. BUCKLE)  THEN 
NFLAG(I) r: 

END  IF 
END  IF 
END  DO 

C  SUMMATION  OF  EACH  TYPE  OF  BREAKAGE 
BURNTaO . 0 
TENSILsO .0 
BUCK=0 . 0 
DO  I3:1,NE 

IF(NFLAG(I) .EQ.l)  BURNT-BURNT+1 
IF (NFLAG (I) .EQ. 2)  TENSIL=TENSIL+1 


74 


on  noononoooono  o  ooo  on  on 


IF (NFLAG ( I ) . EQ . 3 )  BUCK=BUCK+1 
END  DO 

TBROKE =BDRNT+TENS IL+BUCK 
WRITE ( 1 1 , 5 ) BURNT* 1 0 0 /NE 
5  FORMAT ( IX, ' PERCENTAGE  BURNT  FIBERS ' , F5 . 1 ) 

WRITE ( 1 1 , 4 ) TENS IL* 1 0 0 /NE 

4  FORMAT { IX, 'PERCENTAGE  BROKEN  IN  TENSILE  MODE',F5.1) 
WRITE { 1 1 , 2 ) BUCK* 1 0 0 /NE 

2  FORMAT ( IX , ' PERCENTAGE  OF  BUCKLED  FIBERS ' , F5 . 1 ) 

WRITE (11, 1)TBROKE*100/NE 

1  FORMAT (IX, 'TOTAL  PERCENTAGE  OF  FAILED  FIBERS ', F5 . 1 ) 
END 


SUBROUTINE  INTERP (XI , X2 , X3 , Y1 , Y2 , Y3 , XX , YY) 

THIS  SUBROUTINE  TAKES  3  X  VALUES,  3  Y  VALUES.  THEN  USING  A  GIVE  XX 
VALUE  PERFORMS  QUADRATIC  INTERPOLATION  TO  DETERMINE  THE  YY  VALUE. 
IMPLICIT  REAL*8 (A-H,0-Z) 

REAL*8  A(3,3) ,B(3) ,C(3) 

A(l, 1) =1 .DO 
A(2, 1) =1 .DO 
A(3,1)=1.D0 
A(l,2) =X1 
A(2,2) =X2 
A(3,2) =X3 
A(1,3)=X1**2. 

A(2,3) =X2**2 . 

A(3,3) =X3**2 . 

B(l)  =Y1 
B(2)  =Y2 
B(3) =Y3 

CALLING  IMSL  LIBRARY  SUBROUTINE  TO  SOLVE  A(3 , 3) *C (3) =B (3) 

TO  DETERMINE  THE  QUADRATIC  COEFFICIENTS 
CALL  DLSARG(3,A,3,B,1,C) 

YY=C  (1)  +C  (2)  *XX+C  (3)  *XX**2  . 

END 

SUBROUTINE  MOD2FL (NN, NSIZE, Ll , L2 , A, B, S , C) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  S(NN,NN),  C (NN) 
************************** 

APPLY  TYPE  2  CONSTRAINTS  TO  FULL  EQUATIONS 
METHOD  AND  SUBROUTINE  BY:  J.  ED  AKIN 


NN  =  THE  TOTAL  SIZE  OF  THE  PASSED 
NSIZE  =  WORKING  SIZE  OF  ARRAYS 
S  »  FULL  SQUARE  MATRIX 
C  =  COLUMN  MATRIX 

L1,L2  =  CONSTRAINED  DEGREES  OF  FREEDOM  SUBCRIPTS  OF  D  MATIRX 
A,B  =  CONSTRAINT  COEFFICIENTS 

D(L1)  +  A*D(L2)  =  B  ,  S*D  *  C 
INITIAL  CALCULATIONS 
S12  =  S(L1,L2) 

Cll  =  C(L1) 

SUBTRACT  A*COLUMN  Ll  FROM  COLUMN  L2 

SUBTRACT  B*COLUMN  Ll  FROM  C 

DO  10  I  =  1, NSIZE 

S(I,L2)  =  S(I,L2)  -  A*S(I,L1) 

S(L2,I)  =  S(I,L2) 

C(I)  =  C(I)  -  B*S(I,L1) 
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nn  n  n  nnno 


S(I,L1)  =  0.0 
S(L1,I)  =  0.0 
10  CONTINUE 

S(L1,L1)  =  1.0 

S(L2,L2)  =  S(L2,L2)  -  A*S12  +  A*A 
S(L2,L1)  =  A 
S(L1,L2)  =  A 
C(L1)  =  B 

C(L2)  =  C(L2)  -  A*C11  +  A*B 

RETURN 

END 


SUBROUTINE  READER (N1 , NNP , TIME , HEIGHT, YIN, TEMPIN, DIAMIN, NSTOP ) 
THIS  SUBROUTINE  READS  IN  THE  NODAL  POINT  VALDES  FROM  THE 
COMBUSTION  CODE 

IMPLICIT  REAL*8 (A-H,0-Z) 

REAL*8  YIN(Nl) ,TEMPIN(N1) , DIAMIN (Nl) 

READ(12,*,END=10) 

READ(12,*,END=10) 

READ(12,*,END=10) 

READ(12,*,END=10)  TIME 
READ(12,*,END=10) 

READ(12,*,END=10)  HEIGHT 
READ(12,*,END=10) 

CONVERT  HEIGHT  TO  INCHES 
HEIGHT=HEIGHT*12 . 0 
DO  Ial,NNP 

READ(12,*,END=10)  K,YIN(I) ,TEMPIN(I) ,DIAMIN(I) 

CONVERT  INPUT  LOCATION (YIN)  AND  DIAMETER (DIAMIN)  TO  INCHES 
YIN(I) -YIN(I) *12.0 
DIAMIN  (I)  .DIAMIN  (I)  *12 . 0 
END  DO 
GOTO  20 
10  NSTOP.l 
20  END 


SUBROUTINE  CHECK (NNP, NSTOP) 

C  QUICK  CHECK  TO  SEE  IF  THE  NUMBER  OF  INPUT  NODAL  POINTS  IS  ODD 
C  SO  QUADRATIC  INTERPOLATION  CAN  BE  USED.  IF  NOT  STOP  PROGRAM 
MCHECK.MOD (NNP , 2 ) 

IF(MCHECK.EQ.O)  THEN 

WRITE(*,2) 'THE  NUMBER  OF  INPUT  NODAL  POINTS  MOST  BE  AN  ODD' 
WRITB(*,2) 'NUMBER  IN  ORDER  FOR  THIS  PROGRAM  TO  USE  QUADRATIC' 
WRITE (*, 2) 'INTERPOLATION.  CHECK  NNP  IN  INPUT  DATA' 

WRITE (*, 2) 'THIS  PROGRAM  IS  STOPPED' 

?  FORMAT  (IX,  A) 

NSTOP-1 
END  IF 
END 
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APPENDIX  C:  SAMPLE  OUTPUT  DATA  FROM  STRESS  CODE 


TIME=  0.000  SECONDS 

FIBER  NUMBER  FIBER  STRESS  (PSD 


1 

0 .OOOOE+00 

2 

O.OOOOE+00 

3 

O.OOOOE+00 

4 

O.OOOOE+00 

5 

O.OOOOE+00 

6 

O.OOOOE+00 

7 

O.OOOOE+OO 

8 

O.OOOOE+OO 

9 

O.OOOOE+OO 

10 

O.OOOOE+OO 

11 

O.OOOOE+OO 

12 

O.OOOOE+OO 

13 

O.OOOOE+OO 

14 

O.OOOOE+OO 

15 

O.OOOOE+OO 

16 

O.OOOOE+OO 

17 

O.OOOOE+OO 

18 

O.OOOOE+OO 

19 

O.OOOOE+OO 

20 

O.OOOOE+OO 

21 

O.OOOOE+OO 

22 

O.OOOOE+OO 

23 

O.OOOOE+OO 

24 

O.OOOOE+OO 

25 

O.OOOOE+OO 

26 

O.OOOOE+OO 

27 

O.OOOOE+OO 

28 

O.OOOOE+OO 

29 

O.OOOOE+OO 

30 

O.OOOOE+OO 

BURNT  FIBERS 

0.0 

PERCENTAGE  BROKEN  IN  TENSILE  MODE  0.0 
PERCENTAGE  OF  BUCKLED  FIBERS  0.0 
TOTAL  PERCENTAGE  OF  FAILED  FIBERS  0.0 


DIAMETER 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0 .8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 


77 


TIME=  8.810  SECONDS 

FIBER  NUMBER  FIBER  STRESS 


(PSD 


1 

2777. 

2 

2236. 

3 

1791. 

4 

1385. 

5 

1022  . 

6 

705.4 

7 

433.8 

8 

206.0 

9 

18.97 

10 

-131.2 

11 

-249.3 

12 

-340.1 

13 

-408.6 

14 

-459.2 

15 

-495.9 

16 

-522.0 

17 

-540.2 

18 

-552.6 

19 

-561.0 

20 

-566.6 

21 

-570.3 

22 

-572.6 

23 

-574.1 

24 

-575.1 

25 

-575.7 

26 

-576.0 

27 

-576.3 

28 

-576.4 

29 

-576.5 

30 

-576.5 

BURNT  FIBERS 

0.0 

PERCENTAGE  BROKEN  IN  TENSILE  MODE 
PERCENTAGE  OF  BUCKLED  FIBERS  0.0 
TOTAL  PERCENTAGE  OF  FAILED  FIBERS 


0.0 


0.0 


DIAMETER 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0 .8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333B-02 

0.8333B-02 

0.8333E-02 

0.8333E-02 
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TIME=  17.810  SECONDS 


lUMBER 

FIBER  STRESS  (PSD 

DIAMETER 

1 

3524  . 

0.8333E-02 

2 

3022  . 

0.8333E-02 

3 

2597  . 

0 .8333E-02 

4 

2190  . 

0.8333E-02 

5 

1805. 

0.8333E-02 

6 

1443  . 

0.8333E-02 

7 

1107. 

0.8333E-02 

8 

797.8 

0.8333E-02 

9 

515.9 

0.8333E-02 

10 

261.7 

0.8333E-02 

11 

34.97 

0.8333E-02 

12 

-165.1 

0.8333E-02 

13 

-339.8 

0.8333E-02 

14 

-490.6 

0.8333E-02 

IS 

-619.4 

0.8333E-02 

16 

-728.2 

0.8333E-02 

17 

-819.1 

0.8333E-02 

18 

-894.3 

0.8333E-02 

19 

-955.8 

0.8333E-02 

20 

-1006. 

0.8333E-02 

21 

-1046. 

0.8333E-02 

22 

-1077. 

0.8333E-02 

23 

-1102. 

0.8333E-02 

24 

-1122. 

0.8333E-02 

25 

-1136. 

0.8333E-02 

26 

-1148. 

0.8333E-02 

27 

-1156. 

0.8333E-02 

28 

-1161. 

0.8333E-02 

29 

-1165. 

0.8333E-02 

30 

-1167. 

0.8333E-02 

BURNT 

FIBERS  0 . 0 

PERCENTAGE  BROKEN  IN  TENSILE  MODE  0.0 
PERCENTAGE  OF  BUCKLED  FIBERS  13.3 
TOTAL  PERCENTAGE  OF  FAILED  FIBERS  13.3 
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26.560  SECONDS 


MBER 

FIBER  STRESS 

(PSD 

DIAMETER 

1 

1389  . 

0.8333E-02 

2 

1519  . 

0 .8333E-02 

3 

1520. 

0.8333E-02 

4 

1392  . 

0.8333E-02 

5 

1174  . 

0.8333E-02 

6 

902.0 

0.8333E-02 

7 

603.2 

0.8333E-02 

8 

297.3 

0.8333E-02 

9 

-3.231 

0.8333E-02 

10 

-290.8 

0.8333E-02 

11 

-561.4 

0.8333E-02 

12 

-812.8 

0.8333E-02 

13 

-1044  . 

0.8333E-02 

14 

-1255. 

0.8333E-02 

15 

-1445 . 

0.8333E-02 

16 

-1616. 

0.8333E-02 

17 

-1768. 

0.8333E-02 

18 

PREVIOUSLY 

BROKEN, 

NFLAG»3 

19 

PREVIOUSLY 

BROKEN, 

NFLAG=3 

20 

PREVIOUSLY 

BROKEN, 

NFLAG=3 

21 

PREVIOUSLY 

BROKEN, 

NFLAG>3 

22 

PREVIOUSLY 

BROKEN, 

NFLAG=3 

23 

PREVIOUSLY 

BROKEN, 

NFLAGx3 

24 

PREVIOUSLY 

BROKEN, 

NFLAGx3 

25 

PREVIOUSLY 

BROKEN, 

NFLAG-3 

26 

PREVIOUSLY 

BROKEN, 

NFLAGx3 

27 

PREVIOUSLY 

BROKEN, 

NFLAGx3 

28 

PREVIOUSLY 

BROKEN, 

NFLAG>3 

29 

PREVIOUSLY 

BROKEN, 

NFLAGx3 

30 

PREVIOUSLY 

BROKEN, 

NFLAGx3 

PERCENTAGE  BURNT  FIBERS  0.0 
PERCENTAGE  BROKEN  IN  TENSILE  MODE 
PERCENTAGE  OF  BUCKLED  FIBERS  56.7 
TOTAL  PERCENTAGE  OF  FAILED  FIBERS 
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APPENDIX  D:  SAMPLE  INPUT  DATA  TO  COBSSUSTION  CODE 


80. EO  TAMB  AMBIENT  TEMPERATURE 

50.4E0  UINF  DIFF.  PRESS.  ACROSS  MEDIUM  LB/FT2  OR  SURF  VELOCITY 

2116. 8E0  PAMB  AMBIENT  PRESSURE  liB/FT2 

1.4E0  TORT  TORTUOSITY  OF  THE  MEDIUM 

0.86805555S555555507E-03  DINT  LUMP  FIBER  DIAMETER  PER  CELL 

0.868055S55555555507E-03  SINIT  DISTANCE  BETWEEN  LUMPED  FIBERS 

0.416666666666666667E-01  XLENGTH  TOTAL  MEDIUM  THICKNESS 

0 .416666666666666667E-01  XINIT  INITIAL  THICKNESS  OF  THE  MEDIUM 

S3 . 34E0  RGAS  THE  GAS  CONSTANT  FOR  AIR 

86.0E0  GCOND  THERMAL  CONDUCTIVITY  OF  GRAPHITE 

0.231E0  SPHTG  SPECIFIC  HEAT  OF  GRAPHITE 

70.3E0  ROG  DENSITY  OF  GRAPHITE 

0.9E0  EMIS  THE  EMISSIVITY  OF  THE  GRAPHITE 

.91E0  PHI  SHAPE  FACTOR  1.0  FOR  SPHERE  .91  FOR  CYLINDERS 

l.OEO  HLENTH  REF  LENGTH  FOR  SURFACE  HEAT  TRANSFER  (NORMALLY  1) 

14085. EO  ENTHAL  THE  ENTHALPY  OF  (HIEMICAL  REACTION 

0.5E0  ORDNTH  THE  ORDER  OF  THE  CHEMICAL  REACTION 

0.5E0  EXTEMP  THE  EXPONENT  OF  THE  TEMPERATURE  FOR  REACTION  RATE 

0.375E0  STIOCH  STOICHIOMETRIC  RATIO  FOR  THE  REACTION 

2.065E6  RECOEFF  ARRHENIUS  EQUATION  COEFFICIENT 

28836. EO  ENERGY  ACTIVATION  ENERGY  DIVIDED  BY  THE  IDEAL  GAS  CONSTANT 

30000. EO  SQ  THE  HEAT  FLUX  (SIMULATING  THE  FIRE) 

4. DO  T  THE  APPROXIMATE  TIME  INTERVAL  FOR  PRINTING  RESULTS 

O.EO  TSTART  THE  STARTING  TIME  OF  THE  PROBLEM 

200.E0  TSTOP  THE  STOP  TIME  OF  THE  PROBLEM 

40. EO  TQ  THE  HEATER  TURN  OFF  TIME  (FIRE  IS  OUT) 

l.OOD-S  HSTART  THE  INITIAL  TIME  INTERVAL  IN  SUBROUTINE  SDESOL 

l.E-10  HMIN  THE  MINIMUM  TIME  INTERVAL  FOR  SUBROUTINE  SDESOL 

0.500  HMAX  THE  MAXIMUM  TIME  INTERVAL  FOR  SUBROUTINE  SDESOL 

0.0010  RMSEPS  ERROR  CRITERIA  FOR  SUBROUTINE  SDESOL 

2500. EO  TMAX  THE  MAX  TEMP  ON  PRINTER  PLOTS  (Y  AXIS) 

50. EO  TMIN  THE  MIN  TEMP  ON  PRINTER  PLOTS  (Y  AXIS) 

.020E0  CMAX  THE  MAX  OXYGEN  CONC.  ON  PRINTER  PLOTS  (Y  AXIS) 

-.OOIEO  CMIN  THE  MIN  OXYGEN  CONC.  ON  PRINTER  PLOTS  (YAXIS) 

O.EO  SLl  SLOPE  FOR  A  LINEAR  PROFILE  FOR  GRAPHITE  TEMP 

O.EO  SL2  SLOPE  FOR  A  LINEAR  PROFILE  FOR  AIR  TEMP 

O.EO  SL3  SLOPE  FOR  A  LINEAR  PROFILE  FOR  CONCENTRATION 

6  MAXDER  MAX  ORDER  EQUATION  SOLVER  USED  IN  SDESOL 
0  IPRT  PRINT  PARAMETER  FOR  SDESOL 

1  ITIME  MULTIPLIER  OF  PRINT  TIME  STEP  ONCE  PROGRAM  IS  IN  SURF  RECESS 
25  NNP  NUMBER  OF  NODAL  POINTS 

0  lELEM  O^AUTO  1-USER  INPUT (SEE  BELOW)  OF  NODAL  LOCATIONS 

0  IDIAM  O-AUTO  1«USBR  INPUT (SEE  BELOW)  OF  PARTICLE  DIAMETER 

0  ISPACE 

0  lUNIF  OaUINF  IS  DIFFERENTIAL  PRESSURE  1=UNIF  IS  SURFACE  VELOCITY 
0  IREAC  0=AUTO  STOICH.  RATIO  OF  CARBON  OXIDES  1=USER  ENTRY  OF  RATIO 
12  IBC  BOUNDARY  CONDITION  PARAMETER 

1  lie  1.AUTO  LINEAR  TEMPO,  TEMPA,  CONC.  02  0=USER  ENTRY (SEE  BELOW) 

0  IQ  LOCATION  OF  HEATER  OaHEATER  AT  ZERO  laHEATER  AT  XLENGTH 

500  NSTOP  MAX  NUMBER  OF  TIMES  ALLOWED  THROUGH  INTEGRATION  LOOPS 
0  ITRIP  LEAVE  ZERO  UNLESS  RESTARTING  PR(X>RAM  IN  SURF  RECESSION  MODE 
800. EO  --GRAPHITE  TEMP  NODE  1] *THESE  3  ITEMS  ARE  REQUIRED  ONLY 
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800. EO  -  AIR  TEMP  NODE  1 . ] *IF  IIC=1.  THEY  MUST  FOLLOW  "ITRIP"  OR 

.0172E0  --OXYGEN  CONC  NODE  1--] ‘MANUAL  ENTRY  OF  NODE  LOCATIONS  IF  USED 
O.DO  ******THE  FOLLOWING  DATA  IS  THE  X  LOCATION  OF  THE  NODES 
.OOIDO  ‘THESE  LOCATIONS  ARE  ONLY  NEEDED  IF  IELEM=1 

.005D0  ‘REQUIRING  NODAL  LOCATION  INPUT.  IF  IELEM=0 

.OIDO  ‘THE  COMPUTER  CALCULATES  A  LINEAR  SPACING  OF  THE 

.025D0  ‘NODES. 

.05D0  * 

.  IDO  * 

.15D0  *  IF  THIS  MANUAL  ENTRY  IS  USED  IT  MUST 

.2D0  ‘  FOLLOW  INPUT  OF  THE  VARIABLE  "ITRIP" 

.25D0  ‘  AND  HAVE  A  VALUE  FOR  EACH  NODAL  POINT 


.3D0 

.35D0 

.4D0 

.45D0 

.5D0 

.55D0 

.6D0 

.65D0 

.7D0 

.738D0 

.8D0 

.85D0 

.9D0 

.95D0 


l.DO  *******eND  of  nodal  point  LOCATION  INPUT“““““*““ 

800. DO  $$$$$$$$START  OF  MANUAL  ENTRY  OF  GRAPHITE  TEMP  FOR  BACH  NODE 
800. DO  $ 

800. DO  $  $$$$$THIS  INFORMATION  ON  TEMPERATURE  AND  CONCENTRATION 

800. DO  $  IS  ONLY  REQUIRED  IF  IIC«0  REQUIRING  MANUAL  DATA  INPUT 

800. DO  $  FOR  EACH  NODAL  POINT.  THE  $$$$$$  DATA  MUST  FOLLOW  DATA 

800. DO  $  FOR  "ITRIP"  IF  NODAL  LOCATIONS  ARE  NOT  BEING  MANUALLY 

800. DO  $  INPUT  (IF  IBLBM«l)OR  ELSE  THIS  DATA  FOLLOWS  NODAL  POINT 

8 00. DO  $  MANUAL  ENTRY  LOCATIONS 

800. DO  $ 

8 00. DO  $  DATA  MUST  BE  GIVEN  FOR  EACH  NODAL  POINT 


800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO 
800. DO  $ 

800. DO  $ 

800. DO  $$$$$$LAST  MANUAL  ENTRY  OF  GRAPHITE  TEMPERATURE%%%%%%%« 
800. DO  $$$$$$START  MANUAL  ENTRY  OF  AIR  TEMPERATURE  BACH  NODE## 
800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 
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800. DO  $ 

800. DO  $ 

800. DO  $ 

300. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $ 

800. DO  $$$$$LAST  MANUAL  ENTRY  OF  AIR  TEMPERATURE  EACH  NODE 
.0075D0  $$$$START  MANUAL  ENTRY  OF  OXYGEN  CONCENTRATION  EACH  NODE 
.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075DO  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $ 

.0075D0  $$$$$LAST  MANUAL  ENTRY  OF  OXYGEN  CONCENTRATION  EACH  NODE 
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APPENDIX  E:  SAMPLE  OUTPUT  DATA  PROM  COMBUSTION  CODE 


31  NODAL  POINTS 


AT  T1ME= 

0.00  S 

MEDIUM  THICKNESS = 

4166666666667D-01  .4166666666667D-01 
NP  X  TEMP-G 


1  O.OOOOOOOOOOOOOOOOE+00 

2  0.1388888888888889E-02 

3  0 .2777777777777778E-02 

4  0.4166666666666666E-02 

5  0.5SS5555SS5555555E-02 

6  0.6944444444444443E-02 

7  0.8333333333333331E-02 

8  0.9722222222222221E-02 

9  O.llllllllllllllllE-01 

10  0.1250000000000000E-01 

11  0.1388888888888889E-01 

12  0.1527777777777777E-01 

13  0.16£666€€86666€€6E-01 

14  0.180S555555555555E-01 

15  0.1944444444444444E-01 

16  0.2083333333333333E-01 

17  0.2222222222222222E-01 

18  0.2361111111111111B-01 

19  0.2499999999999999E-01 

20  0.2e38888888888888E-01 

21  0,2777777777777777E-01 

22  0.2916666666666666E-01 

23  0.3055555555555555E-01 

24  0 .3194444444444444E-01 

25  0.3333333333333333E-01 

26  0.3472222222222222E-01 

27  0.3611111111111110B-01 

28  0 .3749999999999999B-01 

29  0.3888888888888888B-01 

30  0 .4027777777777777B-01 

31  0.4166666666666666B-01 


0.4000000000000000E-«-03 
0 .4000000000000000E-t-03 
0 .4000000000000000B-t-03 
0.4000000000000000B-t-03 
0 .4000000000000000E-I-03 
0 .4000000000000000E-(-03 
0.4000000000000000E-«-03 
0 .4000000000000000E-t-03 
0 .4000000000000000E-t-03 
0 .4000000000000000E-«-03 
0 .4000000000000000E-«-03 
0.4000000000000000E-t-03 
0 . 4000000000000000E4-03 
0 .4000000000000000E+03 
0 .4000000000000000E+03 
0.4000000000000000B+03 
0 .4000000000000000B+03 
0 .4000000000000000E-(-03 
0 . 4000000000000000E+03 
0 .4000000000000000Ei-03 
0 .4000000000000000E-t-03 
0 .4000000000000000E+03 
0 .4000000000000000B-«-03 
0.4000000000000000E+03 
0 .4000000000000000E-)-03 
0 . 4000000000000000B-t-03 
0.4000000000000000B-<-03 
0 . 40000000000000006-^03 
0 .4000000000000000E-h03 
0 .4000000000000000E-h03 
0 .4000000000000000E-t-03 


DIAMETER 

0 . 69444444444444446-03 
0 .6944444444444444E-03 
0 .6944444444444444E-03 
0 . 6944444444444444E -03 
0 .6944444444444444E-03 
0 .6944444444444444E-03 
0.6944444444444444E-03 
0 .6944444444444444E-03 
0 .6944444444444444E-03 
0. 6944444444444444E-03 
0 .6944444444444444E-03 
0 .6944444444444444E-03 
0.6944444444444444E-03 
0.6944444444444444E-03 
0.6944444444444444E-03 
0 . 6944444444444444E- 03 
0 . 6944444444444444E- 03 
0.6944444444444444E-03 
0.6944444444444444E-03 
0 .6944444444444444E-03 
0 .6944444444444444E-03 
0 . 6944444444444444E- 03 
0.6944444444444444E-03 
0.6944444444444444E-03 
0 . 6944444444444444E- 03 
0 . 6944444444444444E- 03 
0 . 6944444444444444E- 03 
0 .6944444444444444E-03 
0 . 6944444444444444E -03 
0.6944444444444444E-03 
0 . 6944444444444444E- 03 
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AT  TIME= 

8.81  S 

MEDIUM  THICKNESS^ 

4166666666667D-01  .4166666666667D-01 
NP  X  TEMP-G 


1  O.OOOOOOOOOOOOOOOOE+00 

2  0.1388888888888889E-02 

3  0.2777777777777778E-02 

4  0.41666666e6e66666E-02 

5  0.5555555555555555E-02 

6  0 .6944444444444443E-02 

7  0.8333333333333331E-02 

8  0 .9722222222222221E-02 

9  O.llllllllllllllllE-01 

10  0 .1250000000000000E-01 

11  0.1388888888888889E-01 

12  0.1527777777777777E-01 

13  0.1666666666666666E-01 

14  0.1805555555555555E-01 

15  0.1944444444444444E-01 

16  0 .2083333333333333E-01 

17  0.2222222222222222E-01 

18  0.2361111111111111E-01 

19  0.2499999999999999E-01 

20  0.26388a8888888888E-01 

21  0.2777777777777777E-01 

22  0.2916666666666666E-01 

23  0.3055S55S5S5S5S55E-01 

24  0.3194444444444444E-01 

25  0.3333333333333333E-01 

26  0.3472222222222222E-01 

27  0.3611111111111110E-01 

28  0.3749999999999999E-01 

29  0.3888888888888888E-01 

30  0.4027777777777777E-01 

31  0.4166666666666666E-01 


0.6171655414743030E+03 
0 .5781984237740606E+03 
0.5511844753523755E+03 
0.5261308274882458E+03 
0.5035640136879362E-t-03 
0 .4836232366654954E-t-03 
0.4663593795839235E+03 
0 .4517249092494148E4-03 
0 .4395833460910138E-)-03 
0.4297262006976060E-t-03 
0 .4218950356130471E+03 
0 .4158052068913094E-t-03 
0 .4111679694792949E+03 
0 .4077084937514894E+03 
0 .405178S413371006E+03 
0.4033636888937181E-t-03 
0.4020858170994085E+03 
0 .4012020198152163E-I-03 
0 .4006011841982683E-(-03 
0.4001993564763142E-»-03 
0 .39993475S6308410E-t-03 
0.3997630187488629E-I-03 
0.3996530135102078E+03 
0.3995833611737365B+03 
0 .3995396822558629E-(-03 
0 .3995124994579232E+03 
0.3994956957918823E+03 
0 . 3994854173856393E-)-03 
0 . 3994793187448693E+03 
0.3994760668482828E-(-03 
0 .3994750372395511E+03 


DIAMETER 

0 .6944444021540727E-03 
0.6944444306055943E-03 
0 .6944444380346299E-03 
0.6944444413489135E-03 
0 . 6944444428425697E-03 
0.6944444435425437E-03 
0.6944444438859773E-03 
0.6944444440630632E-03 
0.6944444441590080E-03 
0.6944444442134259E-03 
0.6944444442455372E-03 
0.6944444442651040E-03 
0.6944444442773239E-03 
0.6944444442850919E-03 
0 .6944444442900914E-03 
0.6944444442933369E-03 
0.6944444442954588E-03 
0.6944444442968574E-03 
0.6944444442977906E-03 
0.6944444442984252E-03 
0.6944444442988702E-03 
0.6944444442991956E-03 
0.6944444442994469E-03 
0.6944444442996526E-03 
0 . 6944444442998301E-03 
0.6944444442999913E-03 
0.6944444443001414E-03 
0.6944444443002882E-03 
0.6944444443004240E-03 
0.6944444443005773E-03 
0.6944444443006692E-03 
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AT  TIME= 

17.81  S 

MEDIUM  THICKNESS= 

4166666666667D-01  .4166666666667D-01 
NP  X  TEMP-G 


1  O.OOOOOOOOOOOOOOOOE+00 

2  0 .1388888888888889E-02 

3  0.2777777777777778E-02 

4  0 .41666e6666666666B-02 

5  0.555555555555555SE-02 

6  0 .6944444444444443E-02 

7  0 .8333333333333331E-02 

8  0.9722222222222221E-02 

9  O.llllllllllllllllE-01 

10  0.1250000000000000E-01 

11  0.1388888888a88889E-01 

12  0.1527777777777777E-01 

13  0.1666666666666666E-01 

14  0.1805555555555S55E-01 

15  0.1944444444444444E-01 

16  0.2083333333333333E-01 

17  0.2222222222222222B-01 

18  0.2361111111111111E-01 

19  0.2499999999999999E-01 

20  0.2638888888888888E-01 

21  0 .2777777777777777E-01 

22  0.2916666666666666E-01 

23  0.3055555555555555E-01 

24  0.3194444444444444B-01 

25  0.3333333333333333E-01 

26  0.3472222222222222E-01 

27  0.3611111111111110E-01 

28  0.3749999999999999E-01 

29  0.3888888888888888E-01 

30  0.4027777777777777E-01 

31  0.4166666666666666E-01 


0 . 6944483272376198E-t-03 
0 . 6582952483489812E-t-03 
0 . 6330284155199859E-I-03 
0 . 6085590421656422E-)-03 
0.5853068872440880E-t-03 
0 .5633915122869859E+03 
0.5429230303669187E-t-03 
0.5239878470762538E+03 
0 . 5066449982116862E-I-03 
0 .4909237200244683E-(-03 
0.4768226385980613E-I-03 
0 .4643106220267744E-I-03 
0 .4533292103391948E-)>03 
0 .4437963847202381E+03 
0 .4356113l28753018E-t-03 
0 .4286596365129955E-t-03 
0  .4228188585575783E-t-03 
0 .4179634359788205E-t-03 
0 .4139692729896838E-t-03 
0 .4107174182763530E+03 
0 .4080968793508916E-t-03 
0 .4060065623126001E-t-03 
0 .4043564178987307E+03 
0 .4030679223461711E+03 
0 .4020740463757375E+03 
0 .401318a723800529E+03 
0.4007570137582123E+03 
0 .4003529789248906E'f03 
0 .4000805996527675E-»-03 
0.3999226594325125E+03 
0 . 3998706998552636E+03 


DIAMETER 

0.6944440200051797E-03 
0.6944442750792861E-03 
0.6944443547335164E-03 
0. 6944443969 186038E- 03 
0.6944444188805062E-03 
0 .6944444303859835E-03 
0 .6944444364794849E-03 
0.6944444397590113E-03 
0.6944444415612797E-03 
0.6944444425767153E-03 
0.6944444431650298E-03 
0.6944444435161306E-03 
0.6944444437320398E-03 
0.6944444438687254E-03 
0.6944444439576309E-03 
0.6944444440168827E-03 
0.6944444440572177E-03 
0.6944444440851723E-03 
0.6944444441048360B-03 
0.6944444441188342E-03 
0.6944444441288955E-03 
0.6944444441361831E-03 
0.6944444441414949B-03 
0.6944444441453872E-03 
0.6944444441482521B-03 
0.6944444441503701B-03 
0.6944444441519371E-03 
0.6944444441530992E-03 
0.6944444441539341E-03 
0.6944444441545535E-03 
0.6944444441548335E-03 
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AT  TIME= 

26.56  S 

MEDIUM  THICKNESS= 

4166666666667D-01  . 4166666666667D-01 
NP  X  TEMP-G 


1  0  .OOOOOOOOOOOOOOOOE+00 

2  0.1388888888888889E-02 

3  0 . 2777777777777778E-02 

4  0 .4166666666666666E-02 

5  0 .5555555555555555E-02 

6  0 . 6944444444444443E-02 

7  0.8333333333333331E-02 

8  0.9722222222222221E-02 

9  0 . llllllllllllllllE-01 

10  0.1250000000000000E-01 

11  0.1388888888888889E-01 

12  0.1527777777777777E-01 

13  0 .1666666666666666E-01 

14  0.ia05555555555555E-01 

15  0 . 1944444444444444E-01 

16  0 .2083333333333333E-01 

17  0.2222222222222222E-01 

18  0.2361111111111111E-01 

19  0.2499999999999999E-01 

20  0.2638888888888888E-01 

21  0 .2777777777777777E-01 

22  0.2916666666666666E-01 

23  0.30555555S5555555E-01 

24  0.3194444444444444E-01 

25  0.3333333333333333E-01 

26  0.3472222222222222E-01 

27  0.3611111111111110E-01 

28  0 .3749999999999999E-01 

29  0 .3888888888888888E-01 

30  0 .4027777777777777E-01 

31  0 .4166666666666666E-01 


0 .6306595038438272E+03 
0 .6399241320641830E+03 
0.6447192511167083E+03 
0.6406985836760575E+03 
0 .6302556153620338E+03 
0.6156188718308147E+03 
0 .5986698701861527E+03 
0 .5807979040862021E+03 
0 .5629290453016817E+03 
0 .5456325613672719E+03 
0 .5292359938247957E+03 
0 .5139161110240054E+03 
0 .4997594668844450E+03 
0 .4867985002414993E+03 
0.4750316490868472E+03 
0 .4644343493352758E-t-03 
0.4549653630104841E+03 
0 .4465709326513697E+03 
0.4391880102692778E+03 
0 .4327471084903999E+03 
0 .4271749671005763E+03 
0.4223970689333552E+03 
0 .4183399&16653412E-t-03 
0 .4149334943623446E+03 
0.412112S315457292E+03 
0 .409818a488450426S-f03 
0.4080025338059548E-»-03 
0.4066233572360078E+03 
0.4056520065572266E+03 
0.4050713541567405E+03 
0 .4048773671370597E+03 


DIAMETER 

0 .6944429267898032E-03 
0 .6944437360540249E-03 
0 .694444018104  0645E-03 
0 .6944441916810881E-03 
0 .6944442957791966E-03 
0 .69444435736132 18E- 03 
0 .6944443933015021E-03 
0 .6944444141127005E-03 
0 .6944444261547878E-03 
0.6944444331631678E-03 
0 .6944444372871645E-03 
0 .6944444397507130E-03 
0.69444444 12492834E- 03 
0 . 694444442179593 IE -03 
0 .6944444427698287E-03 
0. 694444443152783 IE -03 
0 .6944444434068504E-03 
0 . 6944444435790780E-03 
0 .6944444436982122E-03 
0 .6944444437821565E-03 
0 .6944444438422836E-03 
0.6944444438859626E-03 
0.6944444439180605E-03 
0.6944444439418518E-03 
0.6944444439595736E-03 
0.6944444439727761E-03 
0.6944444439825334E-03 
0.6944444439895999E-03 
0.6944444439944511E-03 
0 .6944444439974 719E- 03 
0 .6944444439985787E-03 
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